0002  AD  -  Final  Technical  Report 


CONTROL  OF  FORMATION  FLYING  SATELLITES 


For  the  Period  May  13, 1999  -  May  12,2000 


by 


Dr.  Peter  M.  Bainum,  Principal  Investigator 
Avaine  Strong,  Research  Asst,  and  Zhaozhi  Tan,  Sr.  Research  Assoc. 

Howard  University 

Department  of  Mechanical  Engineering 
Washington,  D.C.  20059 

E-mail  pbainum@fac.howard.edu 


This  material  is  based  upon  work  supported  by  the 
Defense  Advanced  Research  Projects  Agency 
Defense  Sciences  Office 
DARPA  Order  No.  H899/00 
Control  of  Formation  Flying  Satellites 
Issued  by  DARPA/CMD  under  Contract  #MDA  972-99C-0020 


Any  opinions,  findings  and  conclusions  or  recommendations  expressed 
in  this  material  are  those  of  the  author(s)  and  should  not  be  interpreted 
as  representing  the  official  policies,  either  expressly  or  implied,  of  the 
Defense  Advanced  Research  Projects  Agency  or  the  U.S.  Government. 


DISTRIBUTION  STATEMENT  A 

Approved  for  Public  Release 

Distribution  Unlimited 

VK>  jm o  quality  msm 


20000601  043 


DISTRIBUTION  STATEMENT  AUTHORIZATION  RECORD 


Title:  f )3Q7\fr£i  0/) 


f/ y/ncj  ife/g.  /  / /V'es 


Authorizing  Official:^  fpA  ^h^  - - 

Agency:  ^D_ S-Rj!l£. _  Ph.  No.(/7Q3]_  &  f_ 

j  ~j  Internet  Document:  URL:  _ _ _ 

(DTIC-OCA  Use  Only) 


Distribution  Statement:  (Authorized  by  the  source  above.) 


A:  Approved  for  public  release,  distribution  unlimited. 


[  |  B:  U.  S.  Government  agencies  only.  (Fill  in  reason  and  date  applied).  Other 

requests  shall  be  referred  to  (Insert  controlling  office). 

[  j  C:  U.  S.  Government  agencies  and  their  contractors.  (Fill  in  reason  and  date 
applied).  Other  requests  shall  be  referred  to  (Insert  controlling  office). 

1  1  D:  DoD  and  DoD  contractors  only.  (Fill  in  reason  and  date  applied).  Other 

requests  shall  be  referred  to  (Insert  controlling  office). 

□  E:  DoD  components  only.  (Fill  in  reason  and  date  applied).  Other  requests 
shall  be  referred  to  (Insert  controlling  office). 

|  |  F:  Further  dissemination  only  as  directed  by  (Insert  controlling  DoD  office 

and  date),  or  higher  authority. 

1  1  X:  U.  S.  Government  agencies  and  private  individuals  or  enterprises  eligible 

to  obtain  export-controlled  technical  data  in  accordance  with  DoD  Directive 
5230.25. 


NOTES: 


ar,  He,  /■/> _  /s_  rr  utf 

DTIC  Point  of  Contact  Date 


REPORT  DOCUMENTATION  PAGE 


rumMppav'n 

OMB  No.  07044)188 


1.  ASCNCr  USE  ONLY  (Lu»bknk) 


4  TITLE  AND  SUSHTtf 


t  REMRTMTE 

May  26,  2000 


Final  Technical  Report 

Control  of  Formation  Flying  Satellites 


X  REPORT  TYPE  AND  OATES  COVERED  Final 
.May  13,  1999  -  May  12,  2000 


ft  FUNDING  NUMBERS 

C#MDA  972-99-C-0020 


X  AUTHOfyS) 

Peter  M.  Bainum  ,  Avaine  Strong,  Zhaozhi  Tan 


7.  PERFORMING  ORGANIZATION  KAMELS)  AND  ADDRES${ES) 

Howard  University 

Dept,  of  Mechanical  Engineering 

Washington,  D.C.  20059 


1  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

0002  AD 


X  SPOHSORIMG/llOWTORMfi  AGENCY  NAME{3)  AND  A0DRES$(ES) 

Defense  Advanced  Research  Projects  Agency 
Defense  Sciences  Office 
3701  N.  Fairfax  Drive 
Arlington,  VA  22203-131 


11.  SUPPLEMENTARY  NOTES 


1ft  SPONSORMGMONITORtNG 
AGENCY  REPORT  NUMBER 


12ft  OtSTWBUnOWAVJUUBJUTY  STATEMENT 

DARPA/DSO  -  Dr.  Ephrahim  Garcia 
DARPA/ASBD  Library 

Defense  Technical  Information  Center:  DTIC-BCS 


12ft  DISTRIBUTION  CODE 


13,  abstract  (UgxfrKjm &o wants)  This  report  summarizes  the  work  completed  during  the  reporting 

period  cited  above.  A  technique  for  maintaining  nominal  separation  distance  between 
adjacent  satellites  in  an  elliptically  orbiting  constellaton  is  developed  and  parametric 
trade-off  studies  performed.  This  technique  is  based  on  an  initial  impulsive-type 
correction  of  the  daughter  spacecraft  at  the  first  perigee  in  order  to  shift  the 
direction  of  the  line  of  apsides  by  a  small  angle.  The  technique  works  well  in  Keplerian 
orbits.  In  the  presence  of  perturbations  the  results  are  critically  dependent  on 
the  amplitude  of  the  shift  angle.  Without  subsequent  corrections  near  collision 
situations  can  exist  over  the  long  term.  Additional  feedback  type  of  correctional 
control  is  recommended  to  prevent  unwanted  secular  drifts.  Two  types  of  feedback 
control  strategies  are  considered  here:  an  application  of  the  linear  quadratic  regulator 
(LQR)  theory  based  on  errors  in  position,  and  also  based  on  a  Lyapunov  function  using 
osculating  orbital  elements.  A  preliminary  deployment  strategy  is  introduced  based 
on  near  Hohmann-type  of  transfer  orbits,  and  shows  deployment  can  be  achieved  in 
one  orbit. 


14.  SUBJECT  TERMS  Final  technical  report;  orbital  station  keeping; 

formation  flying;  perturbations;  deployment;  LQR  theory;  Lyapunov 
function 


15.  NUMBER  OP  PAGES 

126 


1ft  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  THIS  REPORT 

unclassified 


NSN  7540*01 -280*5500 


1ft  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

unclassified 


1ft  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

unclass 


2ft  LOMTATTON  OF  ABSTRACT 


Standard  Form  298  (flsv.  2-89) 
Pruett**  ANSI  Sid  Z»1fl 


* 


Y. 

Table  of  Contents 

Page 

Abstract . i 

I.  Introduction  and  Methodology . 1 

II.  Selection  of  Strawman  Configurations  and  Comparison  with  ESA  Auroral 

Cluster  Configuration  and  Proposed  NASA  Multiscale . 14 

III.  Station  Keeping  Maneuvers . 17 

3.1  In  Plane  Station  Keeping  Based  on  Shift  of  Line  of  Apsides . 17 

3.2  Verification  with  Orbital  Dynamics  Simulation  Software . 19 

3.3  Further  Improvement  {consider  thrust  at  apogee  &  perigee} . 26 

3.4  Parametric  Trade  -  Off  Studies  . 38 

3.5  Implementation  -  Propulsive  Requirements  based  on  Lagrange’s 

Planetary  Equations  for  Impulsive  (thrust)  Perturbations . 51 

3.6  Analytical  Solution  to  Station  Keeping  Under  Continuous  Thrust . 55 

3.6.1  Development  of  the  Linear  Quadratic  Regulator  (LQR)  Theory..  56 

3.6.2  Solution  of  the  LQR  based  on  the  Tschauner-  Hempel  Equation  of 

Motion . . 63 

3.6.3  Modification  of  Deployment  Strategy  due  to  the  Earth’s, 

Planetary  and  Lunar  Perturbations . 83 

IV.  The  Implementation  of  Maintaining  Constant  Distance  between 

Satellites  in  Elliptic  Orbits . 88 

4.1  Introduction . 88 

4.2  Review  of  Strategy  for  Maintaining  Distance  in  Elliptic  Orbits . 91 

4.3  Initial  Separation  Deployment . 93 

4.4  Station  Keeping  Maintenance . 98 

4.5  Numerical  Simulations . 105 

4.6  Conclusions . 107 

V.  General  Conclusions  and  Recommendations . 110 

VI.  Implications  for  Further  Research . 112 

References . 115 

Bibliography  of  Papers  Resulting  from  this  Research . 123 


ii 


I.  INTRODUCTION  AND  METHODOLOGY 

The  notion  of  orbital  formation  flying  has  been  known  for  some  time;  recent  budget 
reductions  have  prompted  renewed  interest  in  this  technology.  The  National  Aeronautics  and 
Space  Administration  (NASA)  Cross  Enterprise  Technology  Development  Program 
(CETDP)  is  developing  critical  space  technologies  that  enable  innovative  and  less  costly 
missions[l]  .  Thus  efforts  within  the  Distributed  Spacecraft  thrust  area  (TA)  represent 
technological  developments  which  will  support  NASA’s  ability  to  accomplish  its  task  by 
surpassing  traditional  approaches  to  utilizing  space  and  the  limitation  inherent  to  them,  by 
enhancing  the  ability  of  new  and  emerging  technologies,  and  sciences,  and  further 
accomplishing  the  challenges  of  looking  closer  at  the  world  and  understanding  the  processes 
which  define  in  composition  and  evolution,  looking  deeper  into  the  universe  to  explore 
neighboring  planets,  galaxies  and  to  seek  an  understanding  of  their  origin.  The  trend  is  to 
develop  technology  that  will  allow  a  distributed  network  of  autonomous  spacecraft  to  act 
collaboratively  as  a  single  collective  unit.  An  advantage  of  flying  multiple,  small,  low-cost 
spacecraft  in  formation  is  that  this  provides  correlated  instrument  measurements  formerly 
possible  only  by  flying  many  instruments  on  a  single  large  platform  [2],  Distribution  of 
components  on  a  number  of  satellites  allows  the  advantage  that  a  single  component  failure 
results  in  the  replacement  of  a  small,  cheap  spacecraft  and  not  an  abort  of  the  mission.  Areas 
of  formation  flying  application  in  remote  sensing  include  stereographic 
viewing,interferometry,  conducted  by  DeCou  [3],  and  synthetic  aperture  synthesization  for 

1  Numbers  in  brackets  designate  References  at  end  of  this  report. 
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radar  and  other  applications  conducted  by  Kong,  Miller  and  Sedwick  [4-5].  DeCou 
discussed  the  basic  orbital  configuration  for  interferometry  missions  and  the  thrust 
requirements  for  station  keeping  of  a  two-satellite  formation,  while  Kong,  et  al.,  discussed 
synthesization  of  aperture,  using  orbital  dynamics  which  included  environmental 
perturbations  such  as  non-spherical  Earth  effects,  atmospheric  drag,  solar  radiation  pressure 
and  magnetic  field  interaction.  NASA  has  suggested  at  least  three  formations  of  small 
spacecraft  in  low  earth  orbit  (LEO)  for  scientific  data  collection.  Among  these  are: 

(1)  An  “Auroral  Cluster”  system  of  up  to  four  spacecraft  in  formation  where  the  main 
objective  would  be  to  measure  for  the  first  time  the  curl  of  the  Earth’s  magnetic  field 
vector.  It  is  anticipated  that  this  mission  would  involve  separation  distances  from  a 
few  meters  and  up  to  200  km  over  a  year’s  interval.  For  the  present  study,  500  km 
will  be  considered  the  nominal  separation  distance. 

(2)  A  second  mission  would  involve  two  spacecraft  in  the  same  orbit,  one  following  the 
other,  as  closely  as  5  km,  then  moving  away  to  distances  of  V*  to  Vi  orbit. 
Instrumentation  on-board  would  correlate  distance  measurements  in  the  atmosphere 
over  a  wide  range  of  separation  ranges. 

(3)  A  third  mission  dedicated  to  space  based  interferometry  would  deploy  three 
spacecraft  flying  precise  “zero-drag”  trajectories  with  positional  accuracy  of  the  order 
of  cm,  with  the  objective  of  measuring  gravity  waves  under  zero-drag  conditions.  For 
this  mission  the  interferometry  base  line  set  up  by  the  formation  would  be  of  a  few 
km.  with  the  main  disturbances  attributed  to  the  effects  of  solar  radiation  pressure. 

(4)  In  addition  to  the  three  proposed  Low  Earth  Orbit  (LEO)  missions,  another 
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interesting  concept,  currently  in  the  planning  phase  by  NASA  Goddard  Space  Flight 
Center  (GSFC)  is: 

A  solar  stereo  mission  involving  two  spacecraft  in  heliocentric  orbit,  one  behind  the 
Earth  and  the  second  in  front  of  the  Earth  with  a  separation  distance  of  1 
astronomical  unit  (a.u.)  The  objective  would  be  to  measure  the  solar  corona  mass 
ejection  from  different  vantage  points  within  a  prescribed  base-line  accuracy. 

(5)  In  addition  to  these  proposed  scientific  uses  for  orbital  formation  flying,  clusters  of 
spacecraft  would  be  the  core  of  the  IRIDIUM,  Alcatel,  Globalstar,  Ellipso  and  other 
proposed  orbital  communication  systems  of  the  future.  Depending  on  the  application 
the  number  of  spacecraft  could  vary  from  a  handful  up  to  50  or  1 00s,  with  separation 
distances  from  the  order  of  km  to  thousands  of  km  with  a  wide  range  of  orbital 
altitudes. 

A  search  of  technical  publications,  technical  memoranda,  and  other  reports  pertaining 
to  this  subject  has  been  conducted.  The  majority  of  the  literature  surveyed  on  Constellation 
Station  Keeping,  comes  from  reports  of  completed  research  or  major  phases  of  research 
presented  in  NASA  programs  which  includes  extensive  data  or  theoretical  analyses. 
Furthermore,  these  reports  include  compilations  of  significant  scientific  and  technical  data 
and  information  deemed  to  be  of  continuing  value.  Other  reports  are  mentioned  and 
referenced  throughout  this  document. 

An  early  study  by  Walker  presents  a  systematic  approach  to  the  analysis  of  coverage 
of  the  Earth  by  means  of  circular  orbit  systems  [6].  The  study  ensured  that  every  point  on  the 
Earth’s  surface  can  always  see  at  least  one  satellite  (or  two  satellites  for  double  coverage) 
above  some  minimum  elevation  angle.  It  is  shown  that  five  or  six  (properly  phased) 
geosynchronous  altitude  satellites  can  provide  (continuous  worldwide)  single  coverage,  and 
seven  or  nine  satellites  double  coverage;  and  a  generalized  approach  to  such  coverage 
assessments  is  presented.  Walker  also  recognized  that  low  Earth  orbit  (LEO)  systems  had 
unique  properties  and  advantages,  and  he  proposed  and  analyzed  a  48  satellite  system 
concept  which  became  the  basis  of  the  Loral  Globalstar  system  [7]. 

A  recent  paper  presents  the  Goddard  Space  Flight  Center  (GSFC)  closed-loop  control 
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method  to  fly  in  either  constellations  or  formations  through  the  use  of  an  autonomous  closed 
loop  three-axis  navigation  control  and  innovative  orbit  maintenance  support  [8].  An 
operational  control  method  for  maintenance  of  constellation  formation  flying  is  developed 
in  this  paper.  Examples  were  taken  from  the  Earth  Observing- 1  (EO- 1 )  spacecraft  that  is  part 
of  NASA’s  New  Millennium  Program  (NMP)  and  from  the  proposed  Earth  System  Science 
Program  Office  (ESSPO)  constellations.  Results  can  be  used  to  determine  the 
appropriateness  of  constellations  and  formation  flying  for  a  particular  case  as  well  as  the 
operational  impacts.  After  constellation  and  formation  analysis  was  completed, 
implementation  of  a  maneuver  maintenance  strategy  became  the  driver. 

Brochet  et  al  present  a  new  method  to  solve  the  linear  station  keeping  optimization 
problem  [9].  Several  optimization  models  were  proposed  for  the  problem  of  resetting  drifting 
satellites  to  ensure  desired  coverage  to  include  those  satellites  that  may  fail  during  the  life 
of  the  constellation.  Each  model  consists  in  minimizing  the  total  fuel  consumption  due  to 
maneuvers.  As  in  all  station  keeping  models,  the  objective  is  to  minimize  the  total  fuel 
consumption  associated  with  satellite  maneuvers.  The  minimization  takes  into  account  the 
trajectory  of  each  satellite  and  constraints  on  their  relative  positions.  An  additional  constraint 
is  introduced  to  limit  the  number  of  satellites  that  can  be  simultaneously  controlled.  All 
satellites  cannot  thrust  simultaneously.  This  is  a  Mixed  Integer  Non-Linear  Programming 
(MINLP)  optimization  problem,  containing  Boolean  and  real  variables.  Boolean  variables 
determine  which  satellite  can  be  maneuvered  at  each  step,  and  the  real  variables  correspond 
to  the  thrust  value  of  the  maneuvers.  One  has  to  solve  the  optimization  station  keeping 
problem  or  find  these  values  for  each  cycle.  This  problem  is  split  on  the  basis  of  the 
generalized  Bender’s  decomposition  method  (projection  on  the  Boolean  variable  space)  [10]. 
(1)  The  first  model  is  linear  and  differential  (relative  satellite  positions).  The  subproblem 
(calculation  of  the  impulsive  thrusts)  is  solved  by  a  dual  approach  that  finds  the  solution  in 
a  finite  number  of  steps.  (2)  The  second  model  is  nonlinear  and  non-differential.  It  represents 
the  real  problem  without  simplifications.  The  resolution  of  the  subproblem  is  done  using  a 
direct  search  approach  (Hook  and  Jeeves  algorithm  [1 1])  to  determine  real  variables  in  the 
subproblem.  The  objective  is  to  minimize  the  total  fuel  consumption  of  each  satellite  subject 
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to  (a)  a  coverage  constraint  (distance  between  two  satellites  |  Ad|  must  not  exceed  a  threshold 
O),  and  (b)  three  remaining  operational  constraints. 

In  conclusion[9]  for  the  two  models:  Linear  and  differential  optimization  problem 
and  Nonlinear  and  non-differential  station  keeping  problem,  the  method  presented  is  very 
efficient  since  it  can  provide  the  global  minimum  (optimal  maneuvers  needed  to  ensure  a 
good  coverage)  just  by  resolving  a  system  of  linear  equations  and  can  be  applied  to  every 
station  keeping  problem.  The  second  model  can  solve  more  complicated  cases  (nonlinear 
mixed  variables  problem).  It  is  more  precise  than  the  linear  one.  It  allows  for  various 
constraints  on  position  (i.e.,  relative  or  absolute  positions  of  satellites),  and  various 
definitions  of  the  fuel  consumption.  This  model  can  also  be  used  as  a  first  approach  to 
calculate  optimal  maneuvers  needed  to  replace  a  failed  satellite.  However,  the  global 
problem,  including  the  choice  of  the  spare  satellite,  needs  to  be  stated. 

The  flying  of  spacecraft  in  constellations  and  formations  will  permit  the  accurate 
determination  of  three  -  dimensional  and  time  -  varying  phenomena  and  will  make  it  possible 
to  distinguish  between  spatial  and  temporal  variations.  However,  constellations  and 
formation  flying  impose  additional  complications  on  orbit  selection  and  orbit  maintenance, 
especially  when  each  spacecraft  has  its  own  orbit  or  science  requirements.  Every  object  in 
orbit  follows  approximately  the  second  Keplerian  law  [12].  Therefore,  there  is  a  natural 
tendency  of  separations.  To  maintain  constant  separation  distance  between  pairs  of  spacecraft 
is  not  an  easy  problem  if  the  orbits  are  elliptical.  Hughes  and  Hall  [13]  investigated  an 
optimization  scheme  to  determine  the  best  configuration  for  a  four  spacecraft  formation  in 
a  circular  orbit.  Ulybyshev  [14]  used  the  linear-quadratic  regulator  technique  for  feedback 
control  to  maintain  station  keeping  of  a  circular  orbiting  constellation.  Carpenter  [15] 
suggested  a  distributed  satellite  formation,  modeled  as  an  arbitrary  number  of  fully  connected 
nodes  in  a  network,  for  an  equatorial  circular  Keplerian  orbit,  could  be  controlled  using  a 
decentralized  controller  framework.  Cao,  Modi,  et.  al.,  explored  the  possibility  of  using 
separately  and  a  combination  of  Feedback  Linearization  Technique  (FLT)  and  Linear- 
Quadratic  Regulator  (LQG)  to  study  the  control  theory  of  an  elastic  space  platform-based 
flexible  manipulator  with  four  links,  two  free  to  slew  while  the  other  two  were  permitted  to 
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deploy  [16].  Chao,  Pollard,  and  Janson  [17]  described  a  method  for  determining  cluster 
orbital  elements  and  the  relative  geometry  and  dynamics  of  satellites  under  a  two-body  force 
field  with  the  secular  Earth’s  oblateness  J2  influences(  Because  the  term  J2  will  be  referred 
to  often  in  this  document,  the  reader  is  referred  to  the  derivation,  in  Battin[18].  They  also 
examined  the  disruption  of  the  formation  due  to  natural  perturbations  and  the  feasibility  of 
a  formation-keeping  strategy.  Smith,  Proulx,  et  al.  [19],  investigated  the  use  of  genetic 
algorithms  to  generate  optimal  station  keeping  strategies,  by  constraining  the  orbit  of  an 
Ellipso™  Borealis  satellite  and  developed  the  minimum  fuel  optimal  bum  strategy  by 
minimizing  the  fuel  required  to  maintain  the  orbit  for  a  given  period  of  time.  Sabol,  Bums 
and  McLaughlin[20]  used  the  Draper  Semianalytic  Satellite  Theory  (DSST)  to  study  several 
satellite  formation  flying  designs  and  their  evolution  through  time.  The  DSST  equations  of 
relative  motion  are  used  in  a  manner  similar  to  the  BG14  propagator[21].  The  DSST 
equations  are  derived  for  Keplerian  two  body  dynamics  (like  the  Hill’s  equations,  or  the 
better  known  Clohessy  Wiltshire  (CW)  [22]  equations,  which  are  a  special  form  of  Hill’s 
equations).  BG14  is  a  high  precision  orbit  propagator  simulation  written  by  McDonnell 
Douglas.  The  BG14  propagator  takes  into  account  the  solar  pressure,  oblateness  of  the 
Earth(up  to  40x40)  and  air  density  variation  with  the  Jacchia  70  (J70)  model.  The  DSST 
Semianalytic  Theory  can  also  be  modified  to  include  effects  of  oblateness,  drag,  etc.  To  be 
more  exact,  all  the  references  cited ,  up  to  this  point,  relate  to  these  Hill’s  equations  or  the 
Clohessy-Wiltshire  (C  W)  set  of  equations,  these  have  been  used  primarily  to  analyze  modem 
guidance  and  rendezvous  problems.  These  include  nonlinear  error  analysis,  station  keeping, 
targeting,  surveillance  and  satellite  clustering.  These  equations  have  severe  limitations  and 
drawback.  The  CW  equations  are  of  a  linear  nature  with  constant  coefficients,  and  describe 
the  coasting  terminal  motion  of  the  ferry  vehicle  relative  to  the  target,  when  the  target  is  in 
a  circular  orbit.  The  main  limitation  of  CW  equations  is  that  the  relative  gravitational 
influence  is  expanded  to  only  first  order  terms  in  the  ratio  of  separation  distance  to  orbital 
radius  (p/r)  ending  up  with  a  fast  growing  error  term  as  the  probe  moves  away  from  the 
station  or  target. 

It  may  be  convenient  to  use  the  not  so  well-known  Tschauner  and  Hempel  equations 
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[23].  These  equations  are  sets  of  linear  equations  which  resemble  the  CW  equations  in  their 
derivation  and  describe  the  motion  of  a  spacecraft  near  a  satellite  in  an  arbitrary  elliptic  orbit 
relative  to  a  rotating  orthogonal  coordinate  frame  fixed  in  the  satellite.  Tschauner  and 
Hempel  found  complete  solutions  for  elliptical  orbits  in  terms  of  the  eccentric  anomaly. 
Tschauner  followed  this  article  with  two  others  [24,25].  It  was  a  considerable  task  to  verify 
the  results,  not  solely  due  to  the  fact  that  the  articles  are  written  in  German.  Carter  and  Humi 
[26,27]  applied  these  equations  to  more  general  Keplerian  orbits  and  later  Carter  and 
Brient[28,29,30]  investigated  impulsive  or  thrust  maneuvers  using  these  equations.  A 
complete  search  revealed  only  an  outline  of  the  derivation  of  the  Tschauner  and  Hempel 
equations;  therefore,  a  complete  derivation  is  based  on  Van  der  Ha  and  Mugellesi  [31]. 

It  will  soon  become  apparent,  that  a  multiple  set  of  disciplines  is  needed  to  treat 
the  problem  considered  in  this  report.  These  include: 

•  Physics  -  Orbital  Mechanics 

•  Engineering  -  Linear  and  Nonlinear  Dynamic  Programming  and  Control 
Theory 

•  Mathematics  -  Differential  Equations  and  Mathematical  Modeling 

This  research  consists  of  a  station  keeping  stage  and  a  deployment  stage.  The  totality  is 
summed  up  in  five  chapters.  In  actuality,  the  order  of  events  is  a  deployment  stage  followed 
by  station  keeping.  However,  for  purposes  of  this  research,  it  will  be  assumed  that  the 
desired  configuration  has  been  reached;  therefore,  station  keeping  will  be  analyzed  first,  and 
deployment  will  be  addressed  in  the  later  part  of  this  document.  The  organization  of  this 
report  is  as  follows.  Following  the  Introduction  and  Methodology  (Chapter  One),  is  Chapter 
Two.  Chapter  Two  deals  with  the  Selection  of  the  Strawman  Configurations  and 
Comparison  with  the  ESA  Auroral  Cluster  Configuration  and  the  Proposed  NASA 
Multiscale  System.  Chapter  Three,  Station  Keeping  Maneuvers,  consists  of  two  parts.  The 
first  part  involves,  a  detailed  study  of  station  keeping  of  spacecraft  using  the  BG14 
propagator,  where  instantaneous  impulses  were  conducted  at  perigee  and  apogee  to  cause  the 
in-plane  shift  of  the  line  of  apsides.  Secondly,  the  study  of  station  keeping  in  a  co-planar  non- 
Keplerian  orbit  is  conducted  based  on  the  Tschauner-Hempel  equations  of  motion.  Chapter 
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4  deals  with  the  deployment  and^tation  peeping  of  spacecraft  based  on  a  Lyapunov  function 
control  strategy.  Although  deployment  is  an  important  part  of  formation  flying,  this  subject 
apparently,  until  now,  has  received  little  attention.  A  preliminary  study  by  Badesha,  et  al. 
[32]  involves  the  deployment  of  six  micro-satellites,  one  at  a  time,  from  a  bus.  At  the 
initialization  state,  the  satellites  fly  in  an  along-track  trajectory  separated  by  nominal  spacing. 
The  study  entailed  the  development  of  a  two-body  (bus  and  satellite)  relative  motion 
propagator  based  on  the  Clohessy-Wilshire  equations  with  drag,  from  which  the  relative 
motion  of  the  micro-satellites  is  deduced.  Their  investigation  did  not  result  in  an  optimum 
deployment  strategy.  Badesha  is  hopeful  that  information  gained  will  be  useful  for  future 
studies  involving  n^inimizing  global  fuel  and  cost.  Other  references  pertaining  to 
deployment/Station  Keeping  can  be  found  in  Chapter  4  of  this  document. 

As  a  background,  control  theory  deals  with  the  dynamic  response  of  a  system  to 
commands  or  disturbances.  The  translation  of  physical  design  objectives  into  a  dynamical 
system  gives  rise  to  the  control  problem.  The  key  elements  involved  are: 

•  A  dynamical  system  which  is  to  be  “controlled” 

•  A  desired  output  or  objective  of  the  system 

•  A  set  of  allowable  (or  admissible)  “controls”  (i.e.,  inputs) 

•  A  performance  functional,  (or  cost  function)  which  measures  the 
effectiveness  of  a  given  “control  action” 

The  notion  of  performance  functional,  performance  index  or  cost  function  is  no  stranger  to 
physics.  As  a  brief  review  and  motivation  on  this  topic,  recall  the  statement  from  Hamilton’ s 
principle  for  conservative  systems  in  classical  physics,  “of  all  possible  paths  along  which  a 
dynamical  system  may  move  from  one  point  to  another  within  a  specified  time  interval 
(consistent  with  any  constraints),  the  actual  path  followed  is  that  which  minimizes  the  time 
integral  of  the  difference  between  the  kinetic  and  potential  energies.”  Therefore,  the  cost 
function  or  performance  index  is  defined  as  the  action 

T 

J  =  \L(q,u)dt  (1.1) 


8 


where 

q  =  generalized  coordinate  vector 

(1.2) 

u  -  q  =  generalized  velocities 

(13) 

U(q)  =  potential  energy 

(1.4) 

T(q,u)  =  kinetic  energy 

(1.5) 

L(q,u)=  T(q,u)~  U(q),  the  Lagrangian  of  the  system  (1.6) 
Similarly,  the  general  system  dynamics  are  given  by  the  physics  of  the  problem  and  can 
be  represented  as,  x  =  f(x,u)  (1.7) 


with 

xeR"  «€/?"•  (>•*) 

(superscript  on  f  means  it  can  in  general  be  time- varying) 
while  the  performance  index  ( PI ) 

1  1  * 

J=  —  xTSx+  —  J(xr(£t+  utRu)  (1.9) 

2  2  to 

J  is  defined  over  the  interval  of  interest  [i,N],  where  S  is  an  n  x  n  symmetric  ,  positive 
definite,  constant  matrix,  Q  is  an  n  x  n  symmetric,  positive  definite,  time-varying  matrix  and 
R  is  an  m  x  m  symmetric  time- varying  positive  definite  matrix.  R  is  to  be  fixed  at  the  outset 
by  the  designer  and  the  positive  definite  requirement  is  to  assure  each  element  of  u  is 
bounded.  Although  the  above  equations  are  shown  for  continuous-linear  quadratic  regulator 
systems  they  can  easily  be  set  up  for  a  discrete-  linear  quadratic  regulator[14].  Distinctions 
are  noted  when  dealing  with  non-linear  systems;  as  is  the  case  when  it  is  necessary  to 
construct  a  scheme  that  makes  the  system  follow  (or  track )  a  desired  trajectory  over  some 
time  interval. 

Control  signals  in  physical  systems  are  usually  obtained  from  equipment  which  can 
provide  only  a  limited  amount  of  force  or  energy.  This  necessitates  placing  restrictions  or 
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constraints  upon  the  set  of  controls  which  satisfy  the  constraints  (set  of  admissible  controls). 
For  most  physical  systems,  the  desired  objective  can  be  attained  by  many  admissible  inputs, 
each  of  which  results  in  a  different  response.  Therefore,  the  best  one  needs  to  be  selected. 
This  requires  the  use  of  a  performance  criterion  (functional). 

Linear  programming  is  concerned  with  the  science  of  decision  and  its  application. 
The  concept  of  optimization  is  now  well  rooted  as  a  principle  underlying  the  analysis  of 
many  complex  decision  or  allocation  problems.  It  offers  a  certain  degree  of  philosophical 
elegance  that  is  hard  to  dispute,  and  it  often  offers  an  indispensable  degree  of  operational 
simplicity.  Using  this  optimization  philosophy,  complex  decision  problems  involving  the 
selection  of  values  for  a  number  of  interrelated  variables  will  be  approached  by  focusing 
attention  on  a  single  objective  designed  to  quantify  performance  and  measure  the  quality  of 
the  decision.  The  one  objective  is  maximized  (  or  minimized  )  subject  to  the  constraints  that 
may  limit  the  selection  of  decision  variable  values  [33,34,35].  Stated  mathematically, 
problems  such  as  these  belong  to  the  calculus  of  variations  or  the  Brachistochrone  problem. 

Optimal  Control  Theory  (a  means  of  determining  or  finding  for  a  given  process  the 
control  that  is  best  in  some  sense)  was  developed  as  a  result  of  attempting  to  solve  the 
dynamics  of  classical  systems  which  were  no  longer  linear  time  -  invariant;  that  is,  those 
systems  whose  components  could  no  longer  be  adequately  described  by  ordinary  linear 
differential  equations  with  constant  coefficients.  It  provides  the  means  of  combining  linear 
programming  with  control  theory.  It  is  a  branch  of  applied  mathematics  and  also  of  control 
engineering.  Methods  for  designing  optimal  control  systems  require  sophisticated 
mathematical  tools.  Linear  Quadratic  Regulator  (LQR)  theory  [36]  will  play  major  roles  in 
the  development  of  algorithms  used  in  this  research.  Kluever  and  Tanck  [37]  investigated 
a  feedback  control  law  for  autonomous  station  keeping  maneuvers  based  on  the  discrete 
version  of the  linear  Clohessy-Wiltshire  equations  of  motion  and  the  discrete  time  asymptotic 
LQR.  Tan,  Bainum  and  Strong  [38]  investigated  a  nonlinear  feedback  control  law,  based  on 
a  Lyapunov  function,  for  maintaining  separation  distances  between  several  spacecraft  in  a 
coplanar  elliptical  orbit,  by  applying  this  function  to  the  osculating  orbital  elements  of  the 
daughter  spacecraft.  Typically,  to  solve  a  LQR  problem,  one  begins  by  defining  a 
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performance  index  to  be  optimized  or  the  Hamiltonian  function,  such  as 


H= —(xTQx+ uTRu)+ XT(Ax+ Bu)  (1.10) 


where  the  state  and  costate  equations  are 


dX 

dx 


=  A+  Bu 


-  Qx\  AT X 


(l.ii) 


and  the  stationarity  condition,  provided  u  is  unconstrained  is 

0=  Ru+  BT X  (1.12) 

du 

accordingly  u-  -R~lBTX  (113) 

so  the  optimal  control  sequence  is  determined  if  the  costate  sequence  can  be  found.  The 
constraint  or  system  equation  is  a  recursion  for  the  state  x  that  develops  forward  in  time, 
while  the  costate  equation  is  a  recursion  for  X  that  develops  backward  in  time.  The 
(fictitious)  Lagrange  multiplier  is  thus  a  variable  that  is  determined  by  its  own  dynamical 
equation,  and  it  will  be  seen  that  the  optimal  controller  is  not  causal. 

Deployment  of  spacecraft  is  based  on  maneuvers  from  a  “mother  satellite.”  In  this 
research  two  approaches  are  suggested: 

(1)  Based  on  near  Hohmann  type  transfer  orbits  where  the  impulses  of  the  thrusters  is 
used  to  augment  the  energy  required  from  the  transient  orbit  to  the  final  orbit  and  the 
four  spacecraft  are  deployed  sequentially  from  the  initial  circular  orbit  with  a  small 
time  delay  selected  to  correspond  to  the  desired  shift  in  the  angle  between  the 
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semimajor  axes  of  the  orbits  of  the  adjacent  satellites.  For  the  first  spacecraft  a  single 
impulse  Hohmann  transfer  is  required.  For  the  remaining  three  a  very  small  second 
impulse  would  also  be  required  at  the  second  perigee  position  in  the  final  elliptical 
orbit.  This  method,  described  in  Chapter  4,  is  based  on  the  utilization  of  the  least 
maneuver  energy. 

(2)  Based  on  a  solution  to  the  nonlinear  ( track  a  desired  trajectory  over  time)  two  point 
boundary  value  problem  (TPBVP)  following  Pontryagin’s  Principle[39,40].  The 
TPBVP  is  evident  by  the  state  equation  and  the  adjoint  costate  system,  since  the 
boundary  conditions  required  for  solution  are  the  initial  state  and  the  final  costate 

V-  These  problems  are  in  general  extremely  difficult  to  solve.  Implementation  may 
involve  numerical  techniques  for  solving  TPBVP  such  as  the  quasilinearization  or 
shooting  techniques. 

The  solution  to  two-point  boundary  value  problems  has  been  attempted  in  a  variety 
of  methods,  among  them:  Interpolation  methods,  variational  methods,  method  of  collocation, 
Picard’s  method,  discrete  methods,  finite  different  methods,  quasilinearization,  and  shooting 
methods.  The  ‘shooting  method’  and  quasilinearization  will  be  relied  on  to  solve  these 
equations,  mainly  because  of  the  familiarity  in  these  two  areas[41,42].  Quasilinearization 
will  be  satisfied  with  the  aid  of  recent  advances  in  MATLAB.  A  summary  of  the  two 
techniques  is  listed  below : 

1 .  Quasilinearization.  This  method  is  applicable  only  to  two-point  boundary  value 
problems  for  systems  of  nonlinear  ordinary  differential  equations.  The  original  nonlinear 
problem  is  replaced  by  a  sequence  of  more  easily  solved  linear  problems  whose  solutions 
converge  under  appropriate  conditions  to  the  solution  of  the  original  problem. 

2.  Shooting  methods.  This  method  takes  its  name  from  the  situation  in  the  two-point 
boundary  value  problem  for  a  single  second  order  differential  equation  with  initial  and  final 
values  of  the  solution  prescribed.  The  initial  slope  is  varied  to  give  rise  to  a  set  of  profiles 
which  suggest  the  trajectory  of  a  projectile  “shot”  from  the  initial  point.  That  initial  slope 
is  sought  which  results  in  the  trajectory  “hitting”  the  target,  the  final  value. 

These  methods  may  not  be  so  well-known,  mainly  due  to  perhaps  the  mis-perception 
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that  they  lack  a  certain  amount  of  elegance.  The  reader  can  be  assured  that  these  methods 
contain  all  the  mathematic  rigor  in  development  worthy  of  its  place  in  pure  and  applied 
mathematics  [43,44,45,46]. 

In  this  report  only  the  first  approach  to  the  deployment  problem  will  be  treated.  The 
two  point  boundary  value  problem  approach(TPBVP),  currently  in  progress,  will  be  treated 
in  the  second  year  final  report. 
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II.  SELECTION  OF  BASELINE  (STRAWMAN)  CONFIGURATION 


A  study  was  conducted  of  proposed  NASA  and  ESA  constellation  configurations 
which  would  measure  and  study  upper  atmospheric  phenomena.  The  Auroral  Cluster  (Multi 
scale)  System,  the  Distance  Measurement  System,  the  Orbiting  Interferometer  System,  as 
suggested  by  NASA  for  LEO  mission  together  with  the  Solar  Stereo  System  in  heliocentric 
orbit  were  all  considered  as  possible  baseline  or  “strawman”  configurations  [47].  In  addition 
information  was  also  obtained  from  the  ESA  web  page  on  the  proposed  ESA  Cluster  mission 
with  the  objective  of  determining  physical  processes  involved  in  the  interaction  between  the 
solar  wind  and  the  magnetosphere  by  visiting  key  regions  such  as  the  polar  cusps  and  the 
magnetotail.  Up  to  four  Cluster  spacecraft  would  map  in  three  dimensions  the  plasma 
structure  contained  in  these  regions.  Simultaneous  multi-point  measurements  would  also 
allow  differential  plasma  quantities  to  be  derived  for  the  first  time.  Cluster’s  main  goal  is  to 
study  the  small  scale  plasma  structures  in  space  and  time  in  key  plasma  regions:  solar  wind 
and  bow  shock;  magnetopause;  polar  cusp;  magnetotail;  and  the  auroral  zone.  The 
preliminary  design  of  each  cluster  spacecraft  would  be  based  on  spin  stabilization  at  a 
nominal  rate  of  15  rpm,  with  the  in-orbit  configuration  characterized  by  two  5  m  long 
experiment  radial  booms,  four  50m  long  experiment  wire  booms,  and  two  axial 
telecommunications  antenna  booms,  with  the  spacecraft  diameter  of  2.9  m,  height  of  1 .3  m, 
dry  mass  of  550  kg,  propellant  mass  of  650  kg,  and  payload  mass  of  72  kg  [48]. 

After  reviewing  the  candidate  configurations  for  satellite  constellations,  the  Auroral 
Multiscale  Mission  (AMM)  was  chosen  as  a  baseline  or  “strawman”  model  for  this  research. 
The  reason  for  this  selection  is  that  the  AMM  has  the  objective  of  upper  atmosphere  science, 
and  its  configuration  is  relatively  uncomplicated.  A  brief  description  and  sketch  of  this 
system  is  given  in  Fig.2.1  and  Table  2.1.  It  is  assumed  that  this  relatively  simple  system 
would  not  have  autonomous  navigation  capability.  Initially  for  this  study  the  “strawman” 
configuration  would  be  based  on  four  spacecraft  in  the  same  plane,  principally  along  the  orbit 
track.  The  perigee  altitude  is  selected  to  be  600  km,  the  apogee  altitude  is  selected  to  be  8000 
km.  And  the  nominal  separation  distance  between  two  adjacent  satellites  is  taken  to  be  500 
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km. 
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Figure  2.1 
Auroral  Multiscale 


Single  AMM  Spacecraft/Stowed  Configuration 


Launch  Configuration 
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Table  2.1 

Strawman  Configuration 


i 


Mission 

Multi-point  3-dimensional  data  collection  of  auroral  phenomena 

Orbit 

600x8000  km,  83°  inclination* 

Launch  Vehicle 

Taurus  (21 10  vehicle),  Insertion  Mass  of  465.5  kg,  argp=203.1° 

Spacecraft  Size 

40  inch  diameter  (across  flats) 

Spacecraft  Mass 

90  kg  each 

Science  Payload 

lon/Electron  spectrometer,  UV  Auroral  Imager,  Magnetometer,  Electric 

Fields  (3-axis) 

Position  Knowledge 

GPS,  Knowledge  to  100  meters  (3  sigma) 

Attitude  Knowledge 

0.01°  Star  Tracker  (referenced  to  magnetometer),  star-tracker  implementa¬ 
tion.  Spinning  sun-sensor  /  magnetometer  provides  coarse  attitude. 

Power 

Solar  array  capability:  40  watts  orbit  average  power 

Energy  Storage:  Dual  IPACS  Flywheel  momentum  bias  system  used  for 
both  power  and  attitude  control 

*83  degrees  is  the  proposed  inclination  angle  for  the  AMM.  To  simplify  the  problem 
and  calculations  for  initial  positions/velocities  components,  it  is  easier  to  assume  i = 0,  which 
has  been  done  in  this  report. 
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m..  STATION  KEEPING  MANEUVERS 

3.1  In  Plane  Station  Keeping  Based  on  Shift  of  Line  of  Apsides 


251.32  •  ' 
251.19  * 

251. l|  km. 

Apogee 


V=4.256  km/s 


V=8.771  km/s 

\  I\ 


516.74 
515.66  km. 


Fig.  3.1.1  Variable  Velocities  &  Phase  Distances 


Following  and 
expanding  on  the  article  by 
Tan,  Bainum,  and  Strong 
Perigee  [49],  where  it  was  noted 
that  without  any 
perturbation  or  any  control 
in  a  nominal  Keplerian 
orbit,  given  the  velocity  at 
perigee  and  separation 
distance,  then  one  would 
achieve  near  apogee  the 
separations  shown  in 

Fig.3.1.1.  From  Fig.3.1.1  it  is  seen  that  the  separation  distances  at  perigee  are  more  than 
twice  the  separation  distances  at  apogee.  This  is  a  direct  consequence  of  Kepler’s  law  of 
equal  areas  being  swept  out  in  equal  time  by  the  radius  from  the  Earth  to  the  particular 
spacecraft.  To  maintain  constant  separation  distance  in  such  an  orbit,  one  scenario  would 
be  based  on  the  necessity  to  correct  the  orbit  continuously;  the  tremendous  amount  of  energy 
required  makes  this  strategy  unfeasible. 

A  novel  idea  for  implementation  of  the  separation  distances  for  the  AMM  mission 
was  proposed  in  Ref  [49]  and  is  further  explained  in  this  section.  If  there  are  two  satellites 
in  a  constellation,  the  first  satellite  is  appointed  as  a  reference  satellite,  or  mother  satellite. 
The  semi-major  axis  of  the  orbit  of  the  second  satellite  should  be  shifted  by  a  very  small 
angle  ( something  like  1 0 )  with  respect  to  the  orbit  of  the  mother  satellite  in  order  to  achieve 
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the  approximately  constant  separation  distance  (Fig.3.1.2).  The  propulsive  maneuvers 
required  for  such  a  separation  can  be  calculated  based  on  Lagrange’s  equations  for  impulsive 
thrust  [50],  and  is  presented  in  the  Section  3.5.  If  there  are  two  satellites  flying  in  the  orbits 
shown  as  Fig.3.1.2,  both  co-planar  and,  in  the  counter  clockwise  direction,  at  exactly  the 

same  time  that  they  arrive  at  their 
perigees  and  apogees,  respectively, 
the  distance  between  the  two 
satellites  is  defined  as  the  geometric 
distance,  roughly  the  distance  caused 
by  the  shift  of  the  semi-major  axis  of 
the  orbit;  if  the  satellites  are  flying  in 
exactly  the  same  orbit  (shown  in 
Fig.3.1.1),  the  distance  between  the 
adjacent  satellites  is  defined  as  the 
phase  distance,  i.e  .  the  distance 
accounted  for  the  time  in  which  the 
second  satellite  will  reach  the  mother 
satellite’s  current  position.  From  Fig.3.1.2,  it  is  seen  that  the  geometric  distance  at  the 
perigees  is  smaller  than  the  one  at  the  apogees.  The  strategy  here  is  that  if  the  phase  distance 
is  compensated  by  the  geometric  distance,  then  the  resulting  separation  distance  between  the 
adjacent  satellites  is  maintained  to  be  essentially  constant.  The  details  of  the  calculation  are 
presented  as  follows: 

From  Fig.3.1.2,  if  the  angle  between  the  two  semi-major  axes  is  a,  the  distances 
between  satellites  1  and  2  at  perigee  and  apogee  (the  geometric  distances)  are 

2(R  +  perigee  altitude)Sin(a/2)  and  2(R  +  apogee  altitude)Sin(a/2) 

respectively,  where  R  is  the  radius  of  the  Earth.  Let  the  phase  distances  at  perigee  and  apogee 
be  Pp  and  Pa,  respectively  (Fig  3.1.1),  We  try  to  make  the  resulting  distance  at  apogee  the 
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Orbit  of  Sat  / 

/ 

'-Orbit  of  Safc~2 

Fig.  3.1.2  Angle  Between  Orbits  of  Two  SAT 
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same  as  the  one  at  perigee,  e.g.  500  km.  i.e. 


2(R+perigee  altitude)  sin(—)+Pp  =  500  (3..  1.1) 

2 

2  (R  + apogee  altitude) sin(— )+  Pa  =  500  (3..  1.2) 

2 


From  Fig.3.1.1,  it  is  noticed  that 


Pp  =  517.28 
Pa  251.13 


(3.1.3) 


The  solution  of  equations  (3.1.1M3.1.3)  is  a  =1.341 18979436401  degree.  Pp=336.626  km, 
Pa=163 .426  km.  That  means,  satellite  2  arrives  at  its  own  perigee  (or  apogee)  38.396  seconds 
later  than  satellite  1  does.  Note  the  ratio  defined  as  Pp/Pa  has  a  value  of  2.05980688507. 
The  term  ratio  will  be  used  extensively  in  the  parametric  studies,  section  3.4  ( Table  3.4. 1 
(a)  and  (b)). 

3.2  Verification  With  Orbital  Dynamics  Simulation  Software 

It  is  the  intent  of  this  preliminary  treatment  to  illustrate  the  advantage  of  the  in-plane 
station  keeping  based  on  the  (initial)  shift  of  the  line  of  apsides  between  a  mother  and 
daughter  spacecraft  An  example  of  the  constellation  in  an  orbit  described  above,  with 
separation  distance  of 500  km.,  is  now  given.  The  simulation  is  performed  by  MATLAB  for 
a  little  more  than  half  an  orbit  and  BG14  [21]  for  longer  time  response.  BG14  requires 
precision  orbital  calculations,  based  on  orbital  mechanics,  to  be  used  in  part  as  its 
Initialization  Phase.  Simulation  resuits  are  identical,  as  shown  in  Table  3.2.1  and  Figure 
3.2.1(a),  respectively. 
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Table  3.2.1. 

Time  History  of  Distance  between  Adjacent  Satellites 


T  sec 

D  km 

T  sec 

D  km 

T  sec 

D  km 

T  sec 

D  km 

T  sec 

D  km 

T  sec 

D  km 

* 

500.0 

1152 

482.9 

2323 

477.7 

3494 

488.0 

4646 

497.6 

5760 

499.8 

96 

499.9 

1248 

481.4 

2419 

478.2 

3590 

489.0 

4742 

498.1 

5856 

499.5 

192 

499.3 

1344 

480.2 

2515 

478.9 

3686 

489.9 

4838 

498.5 

5952 

499.3 

288 

498.4 

1440 

479.1 

2630 

479.7 

3782 

490.9 

4934 

498.9 

6048 

498.9 

384 

497.1 

1550 

478.1 

2726 

480.5 

3878 

491.8 

5030 

499.3 

6144 

498.5 

480 

495.6 

1651 

477.5 

2822 

481.4 

3974 

492.6 

5126 

499.5 

6240 

498.1 

576 

493.8 

1747 

477.0 

2918 

482.2 

4070 

493.5 

5222 

499.8 

6336 

497.6 

672 

491.9 

1843 

476.8 

3014 

483.2 

4166 

494.3 

5318 

499.9 

6432 

497.0 

768 

490.0 

1939 

476.7 

3110 

484.1 

4262 

495.0 

5414 

500.0 

6528 

496.4 

864 

488.1 

2035 

476.7 

3206 

485.1 

4358 

495.7 

** 

500.0 

6624 

495.7 

960 

486.2 

2130 

476.9 

3302 

486.1 

4454 

496.4 

5568 

500.0 

6720 

495.0 

1056 

484.5 

2227 

477.2 

3398 

487.0 

4550 

497.0 

5664 

499.9 

6816 

494.3 

T  is  the  time  from  the  perigee,  D  is  the  distance  between  adjacent  satellites,  *  is  the  perigee,  T=0; 
**  is  the  apogee,  T=5490.7  sec.  from  perigee. 

Without  perturbations,  such  as  the  atmosphere  drag,  flatness  of  the  Earth,  the  effect 
of  the  magnetic  fields  of  the  Earth,  the  perturbations  from  the  Moon,  the  Sun  and  the  planets, 
and  so  on,  the  resulting  separation  distance  will  be  kept  essentially  constant  forever;  Fig 
3.2.1(b),  over  a  30  day  duration,  gives  credence  to  this  notion.  The  energy  expended  is  only 
to  be  used  to  compensate  for  these  relatively  small  perturbations.  According  to  this  design, 
significant  control  energy  can  be  saved  while  maintaining  relatively  constant  separation 
distances  between  adjacent  satellites  in  an  elliptical  orbit.  In  addition,  the  control  logic  and 
its  implementation  for  constant  separation  distance  within  a  constellation  are  much  easier 
than  for  the  case  of  time  varying  separation  distance.  If  the  separation  is  varying,  that  is,  the 
separation  distance  is  different  at  each  point  of  the  orbit,  the  measurement  of  the  error  of  the 


20 


separation  distance  requires  the  information  of  the  satellite  location  in  the  orbit.  If  the 
separation  distance  is  constant,  the  measurement  requires  only  a  distance  sensor. 

To  study  the  effect  of  perturbations  due  to  a  non-spherical  Earth  and  the  atmospheric 
drag  the  BG  14  program  was  utilized  with  the  same  initial  conditions  used  in  Table  3.2.1. 
The  effects  of  the  Moon  and  the  Sun  are  also  included  but  are  thought  to  cause  much  smaller 
perturbations.  Figures  3.2.1(a)  through  3.2.2(f)  show  the  results  from  this  numerical 
procedure.  There  is  virtually  no  indistinguishable  difference  in  separation  distance  between 
the  case  with  first  order  oblateness  and  drag  as  compared  with  the  case  of  first  order 
oblateness  without  drag.  This  is  attributed  to  the  relatively  high  perigee  altitude  of  600  km 
used  in  this  example.  This  result  is  different  than  that  shown  in  the  paper  [49],  due  mainly 
to  precision  of  significant  figures.  Fig  3.2.2(a) , which  includes  gravity  degree  4,  drag,  and 
solar  radiation  is  the  same  as  Fig  3.2.2(b)  but  over  a  shorter  time.  Fig  3.2.2(c)  includes  the 
perturbations  of  gravity  degree  4,  drag,  and  solar  radiation,  and  shows  the  response  of  Figs 
3.2.2(a)  and  3.2.2(b)  over  a  longer  period  of  time  with  the  one  radial  thrust-impulse  at  the 
first  perigee;  a  secular  perturbation,  where  the  amplitude  of  oscillation  increases  with  time, 
however  the  mean  value  decreases.  Fig  3.2.2(d)  includes  the  perturbations  of  gravity  degree 
4,  and  solar  radiation.  Fig  3.2.2(e)  includes  the  perturbations  of  gravity  degree  4,  and  drag. 
Fig  3.2.2(f)  is  included  to  show  the  effect  of  no  initial  shifting  of  the  lines  of  apsides  in  a 
perturbed  orbit.  In  the  same  article  by  Badesha  [32],  formation  flying  for  circular  orbits  was 
also  addressed.  It  is  not  the  intent  to  compare  ‘apples  with  oranges,’  but  it  can  be  seen  from 
these  examples,  better  results  and  seemingly  better  control  are  achieved  than  that  of  Badesha. 
In  Badesha’ s  paper,  the  spacing  between  two  spacecraft  oscillate  between  10  and  70  m  with 
the  orbital  period.  The  amplitude  of  the  oscillation  decreases  with  time  but  starts  to  increase 
after  9  days.  After  approximately  21  days,  the  satellites  have  a  potential  to  collide.  It  can 
be  concluded  for  both  the  cases  shown  in  Table  3.2.1  and  Fig.3.2.1,  that  subsequent  station 
keeping  requirements  would  be  far  less  than  those  associated  with  the  original  system  shown 
in  Fig.3.1.1.  As  compared  with  the  Keplerian  results  of  Table  3.2.1  or  Figs  3.2.1(a)  and 
3.2.1(b),  using  the  same  initial  orientation  between  the  lines  of  apsides  of  the  mother  and 
daughter  satellites,  it  is  seen  that  the  separation  distance  at  each  subsequent  apogee  is  less 


21 


than  the  required  500  km. 


Fig.  3.2.1  (a  &  b).  30  Day  Time  History  of  Separation  Distance 

between  Mother  and  Daughter  Satellite  Keplerian  Orbit 
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ReltHve  ctetanca  between  spacecraft  1  and  2 


Plot  of  distance  b«tw««n  spacecraft  14  2 


Fig.  3.2. 2(a  &  b).  The  Distance  between  Adjacent  Satellites 
with  Various  Perturbations  Simulated  by  BG14 
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Fig.  3.2.2(c,d  &  e).  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Various  Perturbations 


2‘ 


Plot  of  spacecraft^  a/0  initial  shift 


Time(days) 


Fig.  3.2.2(f).  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 


Figures  3.2.2(c)  through  3.2.2(e)  show  that  the  Earth’s  oblateness  (J2  )  is  perhaps  the 
dominating  factor  causing  the  secular  perturbation.  In  a  perturbed  orbit,  and  allowing  no 
shifting  of  the  lines  of  apsides  and  the  spacecraft  are  allowed  to  be  separated  as  in  Fig.  3 . 1 . 1 , 
the  secular  perturbation  still  exists(see  Fig.  3.2.2f),  but  the  amplitude  of  the  oscillations  is 
much  greater  than  in  Figs.  3.2.2(c)  through  (e).  Another  interesting  observation  that  can  be 
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especially  recognized  from  Figures3. 2.1(a),  3.2.2(a)  and  3.2.2(b),  is  Kepler’s  second  law  in 
effect  (i.e.,  equal  areas  in  equal  time).  From  these  plots,  perigee  being  the  first  amplitude, 
reveals  a  smaller  area  under  the  curve  than  the  area  under  the  next  amplitude  (apogee).  As 
the  secular  drift  continues,  a  natural  question  would  be,  “what  happens  after  30  days”,  or 
“will  the  satellites  ever  collide  if  not  further  corrected,  and  if  so,  when?”  These  concerns  are 
addressed  in  Section  3.6.3. 

A  series  of  programs  (software)  were  used  to  aid  in  the  execution  of  BG14  and  to 
generate  plots.  One  program  was  used  to  provide  inputs,  based  on  the  desired  configuration, 
as  part  of  the  initialization  phase  for  BG14.  From  the  BG14  input  file  an  output  file  is 
produced.  The  short  routine,  which  follows,  written  in  MATLAB,  is  used  to  make  many  of 
the  plots  contained  in  this  document: 

[trial_dat,  trial_desc]  =  readbgl4(‘  filename  ’); 

plot(trial_dat(:,l),  trial_dat(:,2)); 

xlabel(trial_desc(l,:)); 

ylabel(trial_desc(2,:));  or  alternate  ylabel(‘  appropriate  comments’ ); 

title(‘  desired  caption’ ); 

text(  enter  coordinates,  ’  \  it  { desired  comments }’); 

legend(  ‘a+b\  ‘a  =  sin(x)’,  ‘b  =  cos(x)’);  legend  may  not  be  necessary 

It  should  be  noticed  that  this  program  requires  the  function  routine  (readbgl4)  in  order  to 
read  the  output  file  produced  by  BG14's  input  file. 

33  Further  Improvement  {consider  thrust  at  apogee  and  perigee} 

In  order  to  compensate  for  those  initial  effects,  as  depicted  in  section  3.2,  additional 
impulses  may  be  desired.  From  the  article  by  Schaub,  et.al.[51],  correcting  a  particular 
orbital  element  (a,  e,  i,  12,  o>,  M)  causes  subsequent  errors  in  other  orbit  elements.  Two 
corrections  schemes  will  be  investigated  for  removal  of  the  secular  perturbation.  A  feedback 
control  system  and  a  scheme  made  popular  by  Battin  will  be  investigated  for  possible 
control  sometime  after  the  impulsive  maneuver.  The  correction  scheme  is  now  introduced 
following  the  derivation  of  Battin,  Chap  11,  Sect.  6  [52],  where  Encke’s  method  is  used  to 
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calculate  a  AV,  in  an  attempt  to  further  improve  station  keeping  after  initially  shifting  the  line 
of  apsides .  Encke’s  method  involves  the  planetary  equations  of  motion.  Orbital  elements  are 
derived  from  these  equations  of  motion.  Thus,  correcting  a  single  orbit  element  has  an  effect 
on  the  equations  of  motion,  in  general.  With  the  aid  of  a  transition  matrix,  a  AV  was 
calculated  which  in  turn  was  used  to  calculate  the  impulse  needed  to  make  the  desired 
maneuver  (due  to  the  nature  of  BG14).  The  result,  as  will  be  seen,  of  this  scheme  produced 
no  added  benefit  as  compared  with  the  single  orbit  correction  scheme  as  initially  used .  To 
proceed,  Neustadt[53],  and  also  Stem  and  Potter[54]  noted  that,  at  most,  n  impulses  are 
needed  if  the  state  space  is  n-dimensional.  Thus  at  most  4  impulses  are  needed  for  this 
problem.  Therefore,  a  second  impulse  was  assumed  at  apogee.  To  develop  the  required 
mathematics  and  transition  matrix,  suppose  that  at  some  time,  t,  the  spacecraft  is  found  to 
deviate  from  the  reference  path  by  an  amount  dr(t)  in  position  and  dv{t)  in  velocity.  One 
must  determine  what  the  velocity  deviation  should  be  for  that  particular  position  deviation 
so  that  the  vehicle  will  arrive  at  the  target  point  at  the  predetermined  or  reference  time  (  ta , 
i.e.,  time  at  apogee  in  the  Keplerian  orbit). 


'  dr(t)' 
\dv(t)j 


( 


o 


A 


It  remains  to  determine  the  partitioned  transition  matrix  where  it  is  noted  that  at  any 
time  t  later  than  tx : 
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R'(t,)=  -RT(t) 
where 


G&t,) 


R*(t) 

W(t) 


V(t)J 


^  dx!  dxx 
Kdv!  dxx 


and 


®M) 


(  dr" 
dv) 


®(t„t)  = 


( dxx!  dx 
\dvx!  dx 


dxj  dv^ 
dvx!  dvj 


ref 


dx!  dvx 
dv!  d\x) 


ref 


sothat  Q(t,t1)0(tl,t)=  I  *('iW  V(O=0 

R(tx)  =  0  and  V(tx)  =  T 

R  and  V  matrices  represent  deviations  in  position  and  velocity  from  the  reference  path. 
Also  rref{t)  and  \ ref  (7)  are  the  position  and  velocity  vectors  at  time  t  for  a  vehicle  in  a 

reference  orbit.  Then  ^(f)  =  r'J  v(tJ  &  d v(t)  =  V* (t)d  v(ta  ) 

eliminating  d\(ta)  one  obtains  the  velocity  vector  V  ( whose  gradient  with  respect  to  r 

is  C*)  and  is  that  velocity  required  at  r  to  reach  the  target  point.  If  a  velocity  correction 
A  V(t)  is  to  be  made  at  this  time,  it  can  be  expressed  as 
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a  v(t)  =  £?v(;*)-^v(r)=  c'  (t)Mt' )  -  d  v(r ) 


where  the  superscripts  -  and  +  are  used  to  distinguish  the  velocity  just  prior  to  a  correction 
from  the  velocity  immediately  following  the  correction.  For  these  calculations  to  remain 
valid,  it  is  necessary  to  restrict  the  magnitude  of  the  deviations  from  the  corresponding 
nominal  values.  Alternately,  one  could  target  any  intermediate  point  rT  such  as  the  point 

on  the  planet’s  sphere  of  influence  through  which  the  reference  orbit  passes  at  the  reference 
time  f  .  Then,  if  r  and  v  are  the  position  and  velocity  of  the  vehicle  at  the  time  the 

correction  is  to  be  made,  these  vectors  can  be  extrapolated  to  the  time  tT ,  using  an  orbital 

integration  technique,  such  as  Encke’s  method,  in  order  to  determine  the  point  rT  where 
the  spacecraft  would  be  found  at  the  reference  time  if  no  corrective  action  were  taken.  By 
calculating  the  conic  arc  connecting  the  position  vectors  r  and  r  j  with  a  transfer  time  of 
tT  -  t  (i.e.,  solving  Lambert’s  problem,  Battin[52]  Chapter  7)  the  conic  velocity  vcl  (Non- 

Keplerian  orbit)  at  r  is  determined.  A  second  conic  arc  connecting  the  spacecraft  position 
vector  r  and  the  desired  target  point  tT  produces  the  conic  velocity  vector  vc2  (Keplerian 

orbit).  The  difference  between  the  conic  velocity  and  the  vehicle’s  actual  velocity  is  a  good 
measure  of  the  effect  of  the  disturbing  perturbations.  If  this  velocity  is  corrected  for  the 
effect  of  perturbations,  the  velocity  necessary  to  reach  the  desired  target  from  position  r  is 
obtained.  Therefore  an  excellent  approximation  to  the  required  velocity  correction  is  just  the 
difference  between  the  two  conic  velocities;  specifically, 

A  V  =  vc2-  vcl 
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In  other  words,  the  A  V  needed  for  a  pursuer  spacecraft  (non-Keplerian  orbit)  and  target 
spacecraft  (Keplerian  orbit)  can  be  approximated,  by  noting  the  small  differences  in  velocity 
of  the  two  spacecraft  with  respect  to  small  differences  in  their  respective  orbits. 

This  computation  may  be  repeated  iteratively  to  achieve  the  desired  degree  of 
convergence.  One  would  repeat  the  algorithm(Fig.  3.3.1)  performing  step  2  through  6  as 
often  as  necessary.  At  present  this  is  very  time  consuming  and  laborious.  Obviously  there 
is  a  need  for  automation,  whereby,  the  process  can  function  or  execute  iteratively. 
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Figure  3.3.1  Computation  Algorithm 
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Employing  the  above  algorithm  with  the  following  velocities  and  calculated  values  as  in 
Table  3.3.1. 

Table  3.3.1 

Keplerian  Orbit  Non-Keplerian  Orbit  AV(km/s)  A  V  (km/s)* 


MX) 

.11455918796536 

.14669367598619 

-.032134488021 

.9605749 

73467 

j>d) 

-4.25483469618012 

-4.2641354338482 

.00930073767 

.2780207 

92948 

MX) 

.00000154838281 

.00001488945842 

-1.334107561  X  10-5 

-.01178 

M2) 

-.05991132583111 

-.02700319251021 

-3.21 100658009  XI  O'2 

.9505185 

73882 

m 

-4.25609922462398 

-4.26659407793047 

.0104948533 

.3106674 

72738 

M2) 

.00000166706955 

.00001526738801 

-1 .360031845  X  lO'5 

.0004025 

95104525 

(Apogee  occurs  at  6.35416666667  X  10  2  days,  values  in  table  occurred  at  t  just  prior  to 
apogee  i.e.  6.3194444444  X  10  "2  days )  Both  spacecraft  maneuvered  based  on  an  inertial 
system.  Note:  *  A  V  =  A  V  /|AF| 
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The  bum  times  are  calculated  from  the  information  in  Section  3.5.  In  that  section,  a  force 
impulse  of 3.49397660879  NS  was  found.  That  is,  an  input  of  1  Newton  for  3.49397660879 
seconds.  This  corresponds  to  AV  =.0388219623 199  m/sec.  In  the  case  at  hand,  spacecraft 
no.l  has  a  AVcl  =0334533887605  m/sec  and  spacecraft  no.2,  AVc2  =.0337816289794 
m/sec.  Therefore  AV  =  .0003282402 1 89  m/sec,  using  ratio  and  proportion,  it  is  determined 
that  force  impulses  of  3.01080498844  Nt-Sec  and  3.04034660815  Nt-Sec  are  needed, 
respectively.  The  magnitude  and  direction  for  velocity  components  are  input  via  ‘direction’ 
into  the  BG14  input  file.  These  directions  are  to  be  entered  as  input,  in  unit  vector  form.  It 
is  worth  noting  that  there  are  different  interpretations  in  terminology  as  it  relates  to  the 
direction  component,  in  the  astronautical  community.  Instead  of  the  normal  right  hand 
system  one  may  be  accustomed  to,  BG14  uses  the  right  hand  system  as  depicted  in  Figure 
3.3.2.  ‘X’  is  radial,  ‘Y’  is  along  track,  and  ‘Z’  is  in  the  opposite  direction  (anti  or  nadir 
direction).  In  other  systems,  one  may  find  radial  as  ‘R’  and  the  term  ‘altitude’  used 
interchangeably,  ‘B’  represents  the  tangential  or  along  track  component  and  ‘N’  represents 
the  component  normal  to  an  orbit  plane.  For  the  case  at  hand,  the  ‘  Y’  direction  is  along 
track,  ‘X’  is  radial,  and  ‘Z’  component  input  will  be  opposite  to  that  as  calculated  in  the 
table.  The  inputs  are  entered  in  the  direction  line  of  code,  as  unit  vectors.  Their  order  being: 
along  track  (Y  component), 
crosstrack(Z  component),  and 
radial  (X  component). 


Y 


Y 

z 


X 


Figure  3.3.2 
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520 


Both  spacecraft  maneuvered 


Figure  3.3.3  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 


Utilizing  the  above  method.  Figure  3.3.3,  shows  the  results  with  both  spacecraft 
subject  to  impulsive  maneuvers  at  the  first  apogee.  There  is  no  significant  difference  in  the 
secular  drift  time  history  as  compared  with  the  results  shown  in  Figs.  3.2.2(c)  -  (e)  based  on 
only  the  first  impulse  taken  at  the  first  perigee.  A  second  study,  considered  the  effect  of  a 
second  impulsive  maneuver  of  the  daughter  spacecraft  alone,  at  the  first  apogee,  after  an 
initial  maneuver  at  the  first  perigee.  Then  the  Mother  Spacecraft  and  its  orbit  will  be 
considered  as  the  reference  or  fixed  orbit.  Now  with  a  revised  algorithm,  Table  3.3.2 
becomes  appropriate.  Figure  3.3.4  shows  the  time  response  simulated  by  BG14  over  a  30 
day  period.  In  this  case  a  AV  =  .124057  m/sec  corresponding  to  an  impulse  of  1 1.1651  NS 
and  is  assumed  to  occur  in  the  radial  direction  at  the  first  apogee.  This  represents  about  three 
times  the  impulse  assumed  for  the  first  radial  maneuver  at  the  first  perigee.  By  comparing 
Figure  3.3.4  and  Figures  3.2.2(c,  d  &  e),  it  can  be  seen  that  the  response  of  Figure  3.3.4  has 
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almost  the  same  characteristics  except  for  a  small  amplitude  (“flattening  out”)  of  the 
maximum  and  minimum  peak  amplitudes.  There  is  also  no  significant  difference  between 
the  results  of  Fig.  3.3.3  and  Fig.  3.3.4. 


Table  3.3.2 

Mother’s  Orbit  Daughter’s  Orbit  AV(km/s)  AV(  (km/s) 


*0) 

.16993707421676 

.29387179316892 

-.123934718952 

.99901280915 

8 

*1) 

-4.2633650433302 

4.25785404506646 

-.00551099826 

.04442304706 

49 

i(l) 

.00001498390224 

.00001563415649 

-1.5025425  XI  O'7 

5.24156854 

877  XI  O'6 

(Apogee  for  Mother  is  5490.7  sec  and  5529.096  sec  for  Daughter) 

AV  =  .124057187071  m/sec  or  an  impulse  of  1 1.1651468364  N  S. 
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Daughter  spacecraft  maneuvered  inerilal  frame  of  reference 


Figure  3.3.4  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 

As  a  final  point,  all  the  above  maneuvers  were  performed  based  on  an  inertial  reference 
frame.  The  study  now  continues  with  a  displacement  in  the  relative  motion  in  a  rotating 
LVLH  (local  vertical,  local  horizontal),  Table  3.3.3,  coordinate  system  where  Spacecraft  no.  1 
is  fixed.  An  explanation  of  an  LVLH  coordinate  system,  follows  that  of  Bond  and  Allman 
[55]. 

Table  333 

Chaser  /  Daughter’s  Orbit  AV(km/s)  AV^km/s) 


Ajc 

-497.30555231311 

-.0910816030378 

-.999846745675 

-.00109369196631 

-.0000002003098 

-.00030046483 

Az 

8.70752230559343 

.0015947843 

1.75067174844  XIO2 
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LVLH  Coordinate  System.  (Time  5459.9999997  seconds) 

AV  =  .0910955638269  m/sec  or  animpulseof8. 19860074442  NS.  A  plot  of  this  maneuver 
( not  included  in  document)  supports  the  fact  that  the  physics  is  the  same,  regardless  of  the 
reference  system.  Thus,  the  results  are  no  different  than  those  of  Figures  3.3.3  and  3.3.4. 
However,  nowa  AV  =  0.091 1  m/s  corresponding  to  an  impulse  of 8.1986NS  was  assumed 
to  occur  in  the  radial  direction  at  the  first  apogee. 

Assuming  that  the  Battin’s  scheme  is  successful,  the  secular  drift  will  again  become 
a  concern  if  spacecraft  are  allowed  to  travel  for  some  additional  time  without  a  subsequent 
correction.  In  the  final  analysis,  it  is  clear  that  a  different  type  of  subsequent  maneuver  is 
needed  other  than  continuing  to  shift  lines  of  apsides.  A  feedback  control  will  be  studied 
later  in  Section  3.6  of  this  document. 
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3.4  Parametric  Trade  -  Off  Studies 

Tables  3.4.1(a)  and  (b)  show  a  summary  of  the  MATLAB  input  in  preparation  for 
using  the  BG14  software  in  order  to  conduct  a  parametric  study  over  a  range  of  nominal 
separation  distances  and  phase  difference  ratios  ( or  alternatively  line  of  apsides  orientation 
shift  angles ).  The  phase  distance  ratio  is  defined  as  the  nominal  ratio  of  separation  distance 
at  perigee  to  the  separation  distance  at  apogee  in  the  nominal  uncorrected  Keplerian  orbit. 
In  this  section,  different  than  that  of  Section  3. 1  and  3.2,  distances  between  spacecraft  and 
their  phase  distance  ratio  were  arbitrarily  chosen  as  input.  In  turn,  corresponding  shift  angles 
and  other  initial  conditions  were  calculated.  That  is,  a  priori,  the  ratio  (2.05980966034)  was 
calculated  based  on  the  separation  between  ‘mother  and  daughter’  at  apogee.  In  Table  3.4.1, 
the  initial  conditions  in  cartesian  coordinates  for  position  and  velocity  components  have  been 
calculated  and  are  used  as  inputs  to  the  BG14  orbital  propagator  software.  Note,  in  keeping 
with  the  BG14  coordinate  System,  the  X  component  of  velocity  is  given  in  an  opposite 
direction  than  that  calculated  by  the  program  used  here. 

A  preliminarily  study  using  six  different  scenarios  or  ordered  pairs  (separation,  shift 
angle  in  degrees)  were  chosen  as  follows:  (600,  1.628),  (520,  1.402),  (520,  1.386),  (520, 
1.2609),(500, 1.2124)  and  (100,  .27147).  For  each  case  simulated,  there  are  two  different 
effects : 

( 1 )  The  situation  where  a  4  x  4  non-spherical  Earth’ s  gravitational  field,  aerodynamic 
drag  and  solar  radiation  effects  are  included . 

(2)  The  situation  where  only  the  non-spherical  Earth’ s  gravitational  field  is  included- 
referred  to  as  “no  perturbations” 

It  is  seen  that  there  is  negligible  difference  in  the  response,  for  the  perturbed  and  non- 
perturbed  cases.  For  a  nominal  separation  distance  of  600  km  and  an  initial  line  of  apsides 
shift  angle  of  2.1 0 ,  both  an  increasing  and  decreasing  secular  shift  in  the  mother  -  daughter 
separation  distance  is  such  that  a  near  collision  situation  occurs  approximately  every  16  - 
17  days.  In  addition,  the  maximum  subsequent  separation  distance  exceeds  20,000  km, 
indicating  that  the  two  spacecraft  are  located  on  opposite  sides  of  the  orbit. 

For  combinations  of  separation  distances  and  line  of  apsides  shift  angles  of 
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(520, 2.2 ),  (520, 2.04),  (520, 1 .8)  and  (500, 1 .8),  show  similar  tendencies  to  those  as  the  600 
km  separation,  but  alternating  between  near  collisions  and  maximum  separation  distance  on 
opposite  sides  of  the  orbit.  The  extreme  sensitivity  between  separation  distance  history  and 
the  initial  line  of  apsides  shift  angle  is  particularly  apparent  by  comparison  of  (520,  2.2  ), 
(520,  2.04  );  representing  a  difference  of  0.16  deg  in  shift  angle.  Although  both  scenarios 
are  clearly  unacceptable,  the  frequency  between  the  near  collision  events  and  the  maximum 
separation  distance  epochs  is  noticeably  different. 

The  parametric  selection  of  1 00  km  represents  a  drastic  change  in  nominal  separation 
distance  and  0.27147  deg  in  line  of  apsides  shift  angle.  After  an  initial  secular  decrease  in 
separation  distance  followed  by  near  collision,  there  is  a  secular  increase  in  the  separation 
distance  to  almost  a  maximum  of 4500  km  at  the  end  of  the  30  days.  Simulation  shows  that 
the  100  km  separation  case  has  an  initial  downward  trend,  which  begins  at  about  day  200, 
and  perhaps  a  zero  relative  separation  after  about  800  days. 

Although  the  results  of  this  parametric  study,  up  to  this  point,  are  far  from  conclusive, 
the  results  are  quite  dramatic.  In  order  to  avoid  near  collision  possibility  (followed  by  a 
maximum  separation  extending  to  opposite  sides  of  the  orbit...chaos),  the  success  of  the 
initial  maneuver  is  critically  dependent  on  the  amplitude  of  the  initial  line  of  apsides  shift 
angle  which  results  from  an  initial  radial  impulsive  thrusting  maneuver. 
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Table  3.4.1  (a) 

Trade-off  Studies  Using  BG-14  Software 


S 

1 

8- 

C/3 

1 

Shift  Angle 

xl  km 

^1  km 

x2  km 

y2 

km 

600 

2.1 

1.6289008847 

6964.102386146758 

599.835765352191 

6978.14 

m 

520 

2.2 

1.4022541 

6978.321029613882 

566.8205583353445 

6978.14 

520 

2.1 

1.4117022666 

6967.595400672765 

519.8930767241359 

6978.14 

0.0 

520 

2.04 

1.3862096862 

6967.516326413141 

519.8914670338451 

6978.14 

0.0 

520 

1.8 

1.2609713096 

6967.133025595134 

519.8834931600558 

6978.14 

500 

2.2 

1.3944447015 

6968.501857349625 

499.9070973826904 

6978.14 

mu 

500 

2.1 

1.3574034351 

6968.390758465615 

499.9049432461227 

6978.14 

nu 

500 

2.04 

1.3328914763 

6968.317652302147 

499.9035120363998 

6978.14 

500 

1.8 

1.2124705669 

6967.963280272655 

499.8964237732126 

6978.14 

qi 

490 

2.2 

1.3665544718 

6968.883448342178 

489.912559507882 

6978.14 

m 

490 

2.1 

1.3302541343 

6968.776750739195 

489.9105318712732 

6978.14 

0.0 

490 

2.04 

1.3062324803 

6968.706540682686 

489.909185415297 

6978.14 

0.0 

490 

1.8 

1.1882202775 

6968.366206991123 

489.9025138004483 

6978.14 

250 

2.2 

.69720944502 

6975.730072346111 

249.9883826038613 

6978.14 

£9 

250 

2.1 

.67868981321 

6975.702306271525 

249.9881142267019 

6978.14 

0.0 

250 

2.04 

.66643446716 

6975.68403525283 

249.987937556445 

6978.14 

0.0 

250 

1.8 

.60622679964 

6975.595467936702 

249.9870491241593 

6978.14 

0.0 

100 

2.2 

.27888233266 

6977.754394017853 

99.99925353798552 

6978.14 

100 

2.1 

.27147459208 

6977.749951986059 

99.99922059295476 

6978.14 

Q3| 

100 

2.04 

.26657252458 

6977.747028638967 

99.99923485956863 

6978.14 

B3 

100 

1.8 

.2424897697 

6977.732858668894 

99.99921363543054 

6978.14 

0.0 

z\  =  z2  =  0 
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Table  3.4.1  (b) 

Trade-off  Studies  Using  BG-14  Software 


Separation 

Ratio 

Shift  Angle 

degrees 

xl  km/s 

yl  km/s 

x2 

k 

m 

Is 

yl  km/s 

600 

2.1 

1.628901 

.751663795138 

8.72682812807 

0.0 

8.75913966311 

520 

2.2 

1.4022541 

.20954545507477 

8.76555905578956 

8.76806334707253 

520 

2.1 

1.4117022666 

.65196146981314 

8.73757305465386 

0.0 

8.76186262409621 

520 

2.04 

1.3862096862 

.65197229681564 

8.73764604822093 

8.76193622093917 

520 

1.8 

1.2609713096 

.65202457063491 

8.73799991601468 

0.0 

8.7622929974399 

500 

2.2 

1.3944447015 

.62698519749379 

8.73991894922326 

EEI 

8.76237945291506 

500 

2.1 

1.3574034351 

.62699985494731 

8.74002158570968 

0.0 

8.76248287512021 

500 

2.04 

1.3328914763 

.62700948319525 

8.74008912662367 

8.76255093185438 

500 

1.8 

1.2124705669 

.62705597084535 

8.74041655778747 

0.0 

8.76288085016669 

490 

2.2 

1.3665544718 

.61450020639866 

8.74111151932122 

0.0 

8.76268458275625 

490 

2.1 

1.3302541343 

.6145140039868 

8.74121012999996 

0.0 

8.76278391824823 

490 

2.04 

1.3062324803 

.61452306827408 

8.74127502151602 

EEI 

8.76284928566164 

490 

1.8 

1.1882202775 

.61456682905557 

8.74158960791318 

0.0 

8.76316616643501 

250 

2.2 

.69720944502 

.31402199545938 

8,76253789991239 

8.76816287833623 

250 

2.1 

.67868981321 

.31402383464524 

8.76256374962109 

0.0 

8.76818877733018 

250 

2.04 

.66643446716 

.31440250448148 

8.7625807596817 

(jjQj 

8.7682058198194 

250 

1.8 

.60622679964 

.31403087096743 

8.762663217904 

0.0 

8.76828843380099 

100 

2.2 

.27888233266 

.12566819543109 

8.76888348495991 

Ql 

8.76978392368624 

100 

2.1 

.27147459208 

.12566829343061 

8.76888762984303 

0.0 

8.76978806954898 

100 

2.04 

.26657252458 

.12566840308474 

8.76889035650879 

m 

8.76979079750519 

100 

1.8 

.2424897697 

.12566882106324 

8.76890357598226 

0.0 

8.76980402161085 

il  =  z2  -  0 ,  x  <  0 

For  completeness,  further  studies  of  the  nominal  500  km  separation  distance  were 


studied  concentrating  on  the  small  range  of  shift  angles  near  1 .3414581297943  degrees. 


The  objective  was  twofold:  (1)  To  determine  if  other  optimal  ratios  exists,  and  (2)  to 
determine  the  level  of  sensitivity.  Table  3.4.2(a,b)  contains  summary  of  the  data  used 
surrounding  the  500  km  analysis. 


42 


Table  3.4.2(a) 

Trade-off  Studies  for  500  km  Separation 


Ratio 

Shift  Angle 

degrees 

xl  km 

yl  km 

x2  km 

y2 

km 

2.2 

1.3944447015 

6968.501857349625 

499.9070973826904 

6978.14 

0.0 

2.1 

1.3574034351 

6968.390758465615 

499.9049432461227 

6978.14 

0.0 

2.04 

1.3328914763 

6968.317652302147 

499.9035120363998 

6978.14 

0.0 

2.05980688507 

* 

1.3414581297943 

6962.3129552193 

499.7697508871 

6978.14 

0.0* 

2.05981 

1.34119108298 

6968.342368768165 

499.9039972175107 

6978.14 

0.0 

2.06 

1.34126967 

6968.342603013903 

499.9040014341221 

6978.14 

0.0 

1.8 

1.2124705669 

6967.963280272655 

499.8964237732126 

6978.14 

0.0 

*  Optimal  accepted  value 


Table  3.4.2(b) 

Trade-off  Studies  for  500  km  Separation 


Ratio 

Shift  Angle 

degrees 

xl  km/s 

y\  km/s 

x2 

km/ 

s 

y2  km/s 

2.2 

1.3944447015 

.62698519749379 

8.73991894922326  | 

0.0 

8.76237945291506 

2.1 

1.3574034351 

.62699985494731 

8.74002158570968 

0.0 

8.76248287512021 

2.04 

1.3328914763 

.62700948319525 

8.74008912662367 

0.0 

8.76255093185438 

2.05980688507 

* 

1.34145812979 

43 

.51916724165269 

8.7527292910018 

0.0 

2.05981 

1.34119108298 

.6270062295741 

8.74006629141183 

0.0 

8.76252792236339 

2.06 

1.34126967 

.62700619826428 

8.74006607505535 

0.0 

8.76252770432113 

1.8 

1.2124705669 

.62705597084535 

8.74041655778747 

0.0 

8.76288085016669 

*  Optimal  accepted  value,  X  <  0 


Plots  from  the  scenario  using  a  ratio  of  2.05981,  which  equates  to  a  change  in  ratio 
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by  .00000311493  from  the  nominal  value,  revealed  similar  results  to  those  already 
mentioned.  That  is,  an  alternation  between  near  collisions  and  maximum  separation  distance 
on  opposite  sides  of  the  orbits.  The  change  in  ratio  equates 

to  1.28861 5989913922  x  1CT6  degrees,  hardly  a  realizable  angle  forpractical  purposes. 
To  add,  the  stated  change  in  ratio  causes  the  following  changes  in  positions  and  velocities: 
Increase  in  x1  by  A  x,  =  .0294 1 354886500&7M 

Increase  in  y,  by  A  y,  =  .13424633041070&/w 

Decrease  in  v^by  A  =. 10783898792 14  Ikm!  sec 

Increase  in  vyl  by  A  vyl  =  .01266299958997&W  /  sec 

Increase  in  vy2by  A  vy2  =. 0075355836366 Ikmt  sec 

There  is  room  for  error.  The  above  values  were  generated  via  a  fourth  degree  polynomial. 
Thus  human  selection  is  needed  to  pick  the  desired  root.  Nevertheless,  this  supports  the  fact 
that  very  small  changes  in  ratios  create  significant  changes  in  the  kinematics  variables. 

To  verify  or  to  access  the  possibility  of  other  optimal  values,  then  one  could 
undertake  the  time  laboring  tasks  of  examining  possible  sequences,  and  thus  determining  the 
range  of  values  for  a  calculated  ratio.  Any  elementary  statistic  text  on  the  subject  of  counting 
techniques  will  state: 

In  a  sequence  of n  events  in  which  the  first  one  has  kl  possibilities,  the  second  event 
has  k2,  the  third  has  k3,  etc.,  the  total  possibilities  of  the  sequence  will  be 

kl*  k2*  ky-kn. 

Therefore,  examining  the  sequence  of  numbers  2.05980688507  and  knowing  that  2.05981 
is  out  of  the  range  for  acceptable  station  keeping,  then  there  exist  at  most 

10t  10*  10*  10t  10*  10*  10=  107  -  1  variations  between  the  sequence  of 
numbers.  The  minus  1  accounts  for  the  one  already  known  value.  A  sufficient  proof  would 
be  for  one  to  begin  with  the  unacceptable  value,  while  performing  necessary  variations  and 
proceeding  toward  the  accepted  value,  or  vise  versa. 
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It  is  now  theorized  that  at  most  one  optimal  shift  angle  exists  for  each  nominal 
separation.  That  is,  there  is  one  unique  ratio  (phase  to  geometry  distances)  per  each  given 
separation,  apogee  and  perigee.  To  determine  the  optimal  angle  for  a  general  separation,  a 
systematic  calculation  as  performed  in  Section  3.1  must  be  executed.  Attempts  to  find  this 
unique  ratio  were  investigated  using  Kepler’s  law  dealing  with  equal  areas.  This  method 
also  proved  to  be  unsuccessful  due  to  an  additional  approximation  of  letting  an  arc  be  equal 
to  that  of  a  straight  line.  To  eliminate  the  additional  approximation,  the  laws  of  momentum 
were  used  to  find  these  angles  and  thus  the  relevant  coordinates.  A  glance  of  the  numbers 
in  Table  3.4.3  will  not  reveal  much  discernment.  However,  a  plot  of  the  relative  separation 
using  BG 1 4  showed  the  difference.  The  relative  separation  of  spacecraft  produced  as  a  result 
of  conservation  of  momentum  is  shown  in  Figures  3.4.1  and  3.4.2  for  100  and  400 
kilometers  separations,  and  is  regarded  as  that  based  on  the  most  accurate  input  coordinates. 
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Table  3.4.3 

Coordinates  of  400  Kilometers  Separation 


coordinates  found  by 

using  Kepler’s  law  of 

equal  area 

a=  1 .073 1572546746  degrees 

xl  (km) 

6978.137 

yl  (km) 

0.0 

xl'(km/s) 

0.0 

yi’  ( 

km/s) 

8.76524 

8946906 

3 

x2  (km) 

6983.33319994463 

y2  (km) 

399.9508485336767 

x2'  (km/s) 

-.50200516804901 

x2’ 

(km/s) 

8.75086 

1666858 

54 

coordinates  found  using 

conservation  of 

momentum 

0=1.0731418694681  degrees 

xl  (km) 

6978.137 

yl  (km) 

0.0 

xl'(km/s) 

0.0 

yl'  ( 

km/s) 

8.77006 

3506186 

98 

x2  (km) 

6968.01020467271 

y2  (km) 

399.87178947311 

x2'  (km/s) 

-0.41542669857296 

y2- 

(km/s) 

8.75896 

9265945 

06 

Using  the  data  for  the  input  coordinates  and  velocities  from  Table  3.4.3  and  assuming  that 
an  initial  impulsive  maneuver  to  shift  the  line  of  apsides  has  just  been  completed  at  perigee. 
Figures  3.4.1  and  3.4.2  show  30  day  time  histories  of  separation  distances  between  mother 
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and  daughter  satellites  for  nominal  separation  distances  of 400km  (Figs  3 .4. 1  a  &  b)  and  1 00 
km  (Figs  3.4.2  a  &  b).  For  each  separation  distance  the  first  plot  (Figs  3.4.1a,  3.4.2a)  shows 
separation  distance  in  a  nominal  Keplerian  orbit  and  verifies  the  concept  of  an  initial 
impulsive  shift  in  the  line  of  apsides;  the  second  plot  (Figs  3.4.1b,  3.4.2b)  shows  the  secular 
drift  attributed  to  oblateness  and  atmospheric  drag.  In  comparison  with  previous  results  for 
the  nominal  500  km  separation,  it  is  observed  in  general  that  the  smaller  the  nominal 
separation  distance,  the  smaller  is  the  amplitude  of  the  oscillation  in  the  amplitude  of  the 
separation  distance  for  the  ideal  Keplerian  orbits.  In  addition  the  amplitudes  of  the 
oscillations  for  the  secular  drifts  caused  mainly  by  oblateness,  is  increased  for  the  larger 
nominal  separation  distances. 


Figure  3.4.1(a)  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Keplerian  Orbit 
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Relative  distance  between  spacecraft  1  and  2 


420 


Non-Keplerian  orbit  for  30 


Time(days) 


Figure  3.4.1(b)  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 
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Relative  distance  between  spacecraft  t  and  2 


Keplerian  orbit  for  30 


Time(days) 


Figure  3.4.2(a)  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Keplerian  Orbit 
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Non-Keolerian  orbit  for 


Figure  3.4.2(b)  30  Day  Time  History  of  Separation  Distance 
between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 


To  remove  the  secular  perturbation,  a  feedback  control  system  will  be  investigated  for 
possible  control  sometime  after  the  impulsive  maneuver.  Feedback  control  will  be  studied 
later  in  Section  3.6  of  this  document. 
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3.5  Implementation  -  Propulsive  Requirements  Based  on  Lagrange’s  Planetary 
Equations  for  Impulsive  (thrust)  Perturbations 


In  order  to  study  the  implementation  of  the  in-plane  station  keeping  strategy 
Lagrange’s  planetary  equations  for  impulsive  perturbations  are  evaluated.  These  equations 
are  taken  from  Ref.  [49]  and  are  available  in  several  texts  devoted  to  celestial  mechanics  and 
are  listed  as  follows: 


dQ 
dt 
a H 
~dt 


nar 


nar 


pVI^ 


N  sinw  cosec  i 


-  N  cos  u 
2 


—  =  ^—\j  1  -e  2[R  sin  0  +B(  cos  0  +cos£)] 
dt  p 


—  -RcosQ+B(l  +-)sin0]+2sin2(,/2)/— 

dt  \ie  p  dt 

or  ^=-^/wI[-Rcos0+fi(l+-)sin0]-cos/— 
dt  pe  p  dt 


ae  sin0+#a  ^  g  ] 

dt  »  vTV  r 

dn  _  3  n  da 

dt  2  a  dt 
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The  left  side  represents  the  time-rate  of  change  for  the  various  orbital  elements  and 
parameters.  The  force  impulse  per  unit  mass  is 

F  =  Rr+Nh+Bhxr 

where  R  represents  the  radial  (altitude)  direction  component; 

N  represents  the  component  normal  to  the  orbit  plane; 

B  represents  the  tangential  (along  track)  component; 
r  -  unit  vector  along  the  outward  radius; 

A 

h  -  unit  vector  in  the  direction  of  the  orbital  angular  momentum  vector; 

A 

h  x  r  is  a  unit  vector  in  the  positive  tangential  direction  (local  horizontal). 

The  following  symbols  represent  the  customary  orbital  elements: 
a:  semi-major  axis; 
e:  eccentricity; 

Q:  longitude  of  the  ascending  node; 
i:  orbit  plane  inclination  angle  with  respect  to  the  equator; 
o>:  argument  of  perigee  (perihelion). 

In  addition  n  represents  the  mean  angular  rate  in  the  orbit  and 

G)  =  Q  +  (O  is  the  longitude  of  perigee;  it  is  measured  first  in  the  ecliptic  plane  to 
the  ascending  node,  and  then  in  the  plane  of  the  orbit  to  the  direction  of  perigee.  E  is  the 
eccentric  anomaly,  0  is  the  true  anomaly,  u=  0  +  co  ,  and  p  is  the  semi-latus  rectum 
p  =  h2/p  =  a  (1  -e2),  and  p  =  Earth’s  gravitational  constant. 

In  Fig  3.5.1,  a  geometric  interpretation  is  provided  to  describe  those  six  commonly 
used  elements  of  the  orbit.  The  Sun  is  at  O,  Ox  points  toward  the  vernal  equinox  and  Oz 
toward  the  north  pole  of  the  ecliptic.  Letting  the  plane  of  the  orbit  cut  the  celestial  sphere 
in  the  great  circle  HP  A,  then  H  is  the  point  where  the  body  in  its  orbit  rises  north  of  the 
ecliptic,  called  the  ascending  node.  The  point  at  which  the  body  crosses  the  ecliptic,  moving 
south,  is  the  descending  node.  The  angle  xOH  is  the  longitude  of  the  ascending  node;  it  is 
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written  as  Cl,  measured  eastward 
around  the  ecliptic.  The  angle 
between  the  plane  of  the  orbit  and 
the  ecliptic  (angle  between 

andf  )  is  called  the  inclination;  it 
is  written  as  i  (  often  times  right 
ascension  and  declination  are  used 
interchangeably  with  celestial 
equator  and  vernal  equinox).  The 
angle  HOP  is  called  the  argument 
of  perihelion;  it  is  written  as  cd. 
For  geocentric  orbits,  the  Earth  is  at  point  O  and  the  ecliptic  plane  is  replaced  by  the  Earth’s 
equatorial  plane. 

Back  to  the  research  at  hand,  one  possible  strategy  that  will  result  in  a  change  of  the 
orientation  of  the  semi-major  axis  without  changing  the  length  of  the  semi-major  axis  nor 
the  eccentricity  is  to  provide  force  impulses  only  in  the  radial  direction  at  the  perigee  and  /or 
apogee  positions,  and  maintaining  B  =  N  =  0.  With  this  strategy  the  mean  orbital  angular  rate 
n  remains  constant  as  well  as  the  values  of  Cl  and  i. 

Using  the  parameters  for  the  strawman  configuration:  a  =  10678.14  km.;  e  = 
0.34650029347059;  n  =  5.69423x10"*  rad/sec.  and  p  =  398603.361  km3  /  sec2  the  value  of 
the  force  impulse  can  be  calculated  from  the  expression  for  dco  /dt.  We  recall  that  R  dt  =  (Rf 
/m)  dt  where  Rf  is  the  actual  thrust  in  N  and  m  is  the  mass  of  each  spacecraft  in  the 
constellation  (m  =  90  kg).  If  it  is  desired  to  cause  a  change  in  the  argument  of  perigee  of  1 
degree  (1/57.3  rad)  then  a  force  impulse  of 3.49397660879  NS  at  perigee  (or  apogee )  would 
be  required.  This  appears  to  be  well  inside  the  performance  limits  for  the  proposed  pulse 
plasma  thrusters  (PPT).  Of  course  other  thrusting  strategies  could  be  employed  at  different 
points  along  the  orbit  and  utilizing  other  components  of  the  thrust  vector.  The  example 
presented  here  is  to  show  capability  and  may  not  represent  an  “optimum”  design  strategy. 

To  convert  impulse  in  terms  of  AV 


Figure  3.5.1 
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F  =  dP/dt 
or 

A  tF=  A/wv+  wA  v 

assuming  mass  is  constant,  then 

A  v=  Fdt/m=  0.0388219623199  m/sec. 
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3.6  Analytical  Solution  to  Station  Keeping  Under  Continuous  Thrust  Via 
Tschauner  -  Hempel  Equation  of  Motion 

In  this  section,  the  relative  motion  of  two  spacecraft  (ml  and  m2)  is  discussed.  Each 
spacecraft  is  orbiting  a  third  (massive)  body  such  as  the  Earth  (m3).  The  general  (linear) 
state  equations  are  replaced  by  the  Tschauner  -  Hempel  equations.  In  the  present  form,  the 
equation 

r+-^-r  =  F  has  a  singularity  at  r  =  0.  F 
r 


represents  the  perturbation  forces,  fl  is  the  reduced  mass  and  t *  denotes  the  radial 


displacement. 


This  equation  has  the  integral 


—  =  const.  , 
x 


or  where  F  =  0  therefore  r  ->  oo  as  r  ->  0  . 


If  the  true  anomaly  of  the  target  vehicle  is  used  as  the  independent  variable  (  as  does 
Tschauner  -  Hempel)  instead  of  time,  then  this  should  eliminate  the  requirements  to 
numerically  solve  Kepler’s  equations  at  each  time  step  as  well  as  accomplishing  the  task  of 
regularizing  the  two  body  problem.  Other  forms  of  regularization  have  been  accomplished 
by  using  a  “fictitious”  time,  the  so-called  Sundmann  transformation[55].  Regularization  is 
defined  as  the  elimination  of  singularities  from  the  differential  equation.  However,  the 
process  of  eliminating  a  singularity  from  a  differential  equation  of  motion  by  the  use  of  one 
or  more  integrals  is  referred  to  as  embedding. 

A  general  formulation  for  the  relative  motion  allowing  for  arbitrary  perturbing  or 


thrust  forces  on  each  of  the  two  spacecraft  is  presented.  F  Has  the  form  F  = 


P- 


dr 


where  the  P  components  are  perturbations  not  generally  derivable  from  a  potential,  such  as 
thrust  and  drag.  V(r,t)  is  the  perturbing  potential  that  can  depend  explicitly  on  time  if 
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required.  Approximate  perturbation  solutions  for  the  linearized  relative  motion  equation 
with  continuous  radial  or  circumferential  forces  acting  on  the  spacecraft  will  be  established. 

A  pseudo  method  for  station  keeping  will  now  be  investigated.  After  a  shift  in  the 
lines  of  apsides  is  made,  the  Tschauner-Hempel  equations  involving  a  continuous  control 
function  will  be  used  to  find  a  control  law  for  the  control  of  thrusting,  using  continuous 
linear  quadratic  regulator  (LQR)  theory.  The  derivation  of  the  T  schauner-Hempel  equations 
is  outlined  in  Ref.[31].  In  obtaining  a  general  analytical  solution,  in  an  earlier  study, 
Tschauner  obtained  a  canonical  transformation  that  simplifies  the  homogeneous  part  of  the 
equation.  The  control  influence  matrix,  B,  ordinarily  becomes  quite  complicated  and  thus 
no  advantage  gained.  However,  the  B  matrix  in  this  case  of  study,  is  also  simplified.  Thus 
a  delightful  added  benefit  gained  by  the  use  of  these  equations,  besides  the  fact  they  are 
already  linearized,  is  that  they  exhibit  a  relative  motion  in  the  local  vertical/local  horizontal 
system. 

3.6.1  Development  of  LQR 

The  most  general  development  of  station  keeping  strategy,  using  feedback,  will  now 
be  developed.  The  adaptation  is  as  outlined  in  Kirk  [56]  and  Lewis  [57],  where  the  interest 
is  not  in  regulating  a  state  near  zero  but  the  state  is  to  follow  a  nonzero  reference.  Therefore, 
the  optimal  control  will  consist  of  a  feedback  and  a  feedforward  command.  Given  this  robust 
development,  restrictions  can  be  imposed  whereby  adaptations  can  be  made  for  almost  any 
situation. 

For  a  fixed  final  state  -  the  desired  system  is  to  maintain  a  constant  state  r(0f)  =  R*, 
or  separation  for  all  0,  and  more  specifically,  some  0  =  0f  (reachability  is  required) 
x(0f)=r(0f)  =  Rn  (3.6.1. 1) 

and  given  the  system  x  (0)  =  A(0)  +  B(0)u  +  d  xeR"  and  ueR”  (3.6. 1 .2) 

d  is  a  small  disturbance;  the  most  general  form  of  quadratic  performance  index  ( PI ) 
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J  —  2^(^/)  Rfj  |  + 

\9/  ' 

-{{  [x(0)-r(0)]  Q{9)[x{e)-r{e)]+uT{0)R{d)u{e)  }d6 
Z  ■ 


2 

S  + 


2 

|M)-K«)||  0+||«Wlf  R  }d$  (3.6. 1.3) 


The  performance  measure  indicates  that  the  state  x(0)  is  to  be  maintained  close  to 

Rn  and  r(0)  without  excessive  expenditure  of  control  effort.  R  is  a  function  of  0 

only  if  it  is  desired  to  vary  the  weighting  on  the  control  effort.  S  and  Q  are  real  symmetric 
positive  semi-definite  matrices  (i.e.,  ^  0),  whereby,  R  is  a  real  symmetric  positive  definite 
matrix  (i.e.,  >  0 ).  This  guarantees  the  existence  of  R'1  for  all  6  €  [^0,^] .  Once  the  PI 

weighting  matrices  Q  and  R  have  been  selected,  the  determination  of  the  optimal  feedback 
gain  K  is  a  formal  procedure  relying  on  the  solution  of  nonlinear  coupled  matrix  equations. 
Therefore,  the  engineering  judgment  in  modem  LQR  design  appears  in  the  selection  of  Q  and 
R.  There  are  some  guidelines  for  this;  all  of  which  will  be  discussed  later.  The  objective  is 
to  find  the  optimal  control  sequence  uradial ,  that  maintains  the  desired 

separation,  x(0^)  =  RN  while  minimizing  the  performance  function. 

Now  as  in  Chapter  1 ,  but  for  a  continuous  system,  the  Hamiltonian  can  be  represented 
as 

H{x{d),u(<0),p{0),d)=  ~|x(0>-  r(0)||20+  |lMI2*+  pT{Ax+  Bu}  (3.6.1.4) 


57 


’  ,  n 
x  -  ~z  =  Ax  +  Bu 

dp 


(3 .6. 1.5) 


3H 

dx 


-  -Qx-  Arp+  Qr 


(3.6. 1.6) 


Instead  of  using  the  stationary  condition 


dH  T 

—  =  0  =  Ru+  Bp 
du 


(3 .6. 1.7) 


u=-R1BTp  (3.6. 1.8) 

This  is  the  minimum  energy  constrained  input  control  expressed  as  a  costate  feedback. 
Substituting  (3.6. 1 .8)  into  the  state  equation  (3.6. 1 .5),  yields  the  state  and  costate  equations 


JC* 

"  A(l,l)  :-BR  lBT(  1,2)' 

JC 

'  d  ' 

p\ 

-2(2,1)  i  -Ar( 2,2) 

_P_ 

T 

lQr\ 

The  solution  to  these  linear  and  0  varying  differential  equations  where  Qr  is  a  forcing 
function  is 


x(9f) 

p(Qf) 


=  <p(0f>O) 


x(0) 

p{9) 


d 


dr 


(3.61.10) 


where  (p  is  the  transition  matrix  of  the  system  (3 .6. 1 .9).  If  <p  is  partitioned,  and  the  integral 
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replaced  by  the  2n  x  1  vector 


(3.6.1.11) 


/,(*) 

A(0) 


then  equation  (3.6.1.10)  can  be  written  as 


x(6f)=  <pn(9f,9)x(9)+(p^(9f9)p(9)+fi(9) 

p(9f)=  <p2l(9f,9)x(.9)+  <p22(9f,9)p(9)+  f2(9)  3.6.1.12(a)&(b) 


Now  seeking  to  investigate  the  boundary  conditions,  Lewis  [57]  offers  an  excellent 
interpretation.  It  is  assumed  that  x(0  0)  is  known.  Rewriting  the  cost  function,  equation 
(3 .6. 1.3)  as 


ef 

J=  F(x(0f,0f))  +  j[l(*,M,0)+  pT(0)(f(x,u,0)-  x')]d0 


=  F(x(0f,9f))+  \[H(x,u,9)- pT(9)x’}l9 


ea 


(3.6.1.13) 


since  H(x,u,0)  =  L(x,u,0)  +  pT  (0)f(x,u,0)  (3.6.1.14) 

Using  Leibniz’s  rule,  the  increment  in  J  as  a  function  of  increments  in  x,  p,  u  and  0  is 


dJ  =  (Fx)Tdx\t/  +  F,d9 |,;  +  ( H - pTx’)de\Sj  -  ( H - pTx')d0\„t 

Of 

+  \\htx0x+  Ht Udu-  pr dx'\{Hp  -  x')rdp^i0 


(3.6.1.15) 
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Subscripts  denote  a  derivative  with  respect  to  that  variable.  To  eliminate  the  variation  of 
x' ,  integrate  by  parts  and  find 


r  dx  'dO  =  -  pTdx 


+  | pT  dxdO 


(3.6.1.16) 


Substitute  the  fact  from  reference[56] 

dx(Of  )  =  dx(Of  )  +  jc'  (Bf  )dOf  (3.6. 1.17) 

into  equation  (3.6.1.16)  then 

dJ  -  (Fx  -  p)T cb^0f  +  (F$  +  H- pTx'+  pTx')dO 
-(H-  pT x'+  pT x')do\eo  +  pTdx0o 

+  } {(Hx  +  p')Tdx+  Hrudu  +  (Hp  -  x')Tdp]dO 


(3.6.1.18) 


Since  RN  is  fixed,  then  dx\ 


Of 


=  0.  Therefore,  the  boundary  conditions  imply 


p(0f)  =  Sx(6f ) -  Sr (0f)  (3.6.1.19) 

Using  equation  (3.6.1.19)  and  substituting  x(0f)  from  3.6.1.12(a)  into  3.6.1.12(b)  then 


s[«>11(0/,<?)*(0)  +  </>nOfO)p(0)  +  AO)]-  Sr(8f) 
=  <p2,O/O)x(.0)+  <p21OfO)pO)+  AO) 


Solving  for  p(0)  yields 
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(3.6.1.21) 


p(0)  =  [«>22(^/,»)-s«>12(tf/,s)]  \s<pu(el  ,6)- <p2t(e,  ,S)^x(0) 
+[<pn(e,,e)-  s<pn(e,,e')]\sue)-  srn  -  /,(«)] 


or  p{0)  s  k(0)x(0)  +  m(0)  (3.6.1.22)  Then 

u(0)=  -[R-\0)BT(0)k(0)x(0)  +  R-'(6)BJ (e)m(6)\  (3.6.1..23) 


u(0)=  -F(0)x(0)  +  v(0)  (2nd  term  in3.6.1.  23  <  0) 

(3.6.1.24) 

The  first  term  (coefficient  of  x  in  equation  3.6. 1 .24),  is  the  feedback  gain  and  the  second  is 
the  command  signal  or  feedforward  gain.  k(0)  and  m(0)  of  equation  (3.6. 1 .23)  are  the  unique 
positive  semi-definite  solutions  of  the  Riccati  equations.  Furthermore,  A(0),  B(0),  R(0){R 
is  not  necessarily  0  dependent},  Q(0)  and  S  must  be  found  or  specified.  By  adjusting  the 
elements’  values,  the  weight  and  relative  importance  of  the  deviation  of  each  of  the  states 
from  their  desired  values  is  understood. 

To  avoid  determining  the  transition  matrix,  differentiate  (3.6.1.22) 

p'{0)  =  k'(0)x{0)  +  k{0)x\0)  +  m’(0)  (3.6.1.25) 


Substituting  from  (3. 6. 1.9)  for  p'(0)  and  x’(0)  and  using  (3.6.1.22)  to  eliminate  p(0). 


then 

[  k :'+Q+  kA+  ATk-kBRTlBTk  ]x 
+  [w'+  A Tm  -  kBR~xBTm  -  Qr  +  kd\  =  0 


(3.6.1.26) 
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(recalling  matrices  are  0  dependent)  also  for  every  x(0)  and  r(0),  then 

-  k'=  kA+  Ark  +  Q-  kBR-'BTk 

r  ,  ,  (3.6.1.27) 

-rri-  \  at  -  kBR  lBT]m-  Qr+  kd 

The  set  of  equations  in  (3.6.1.27)  are  considered  as  differential  Riccati  equations  .  If  the 
output  is  specified  as  a  linear  combination  of  states  or  y{0)  =  Cx(0)  ,  then  the  second 

term  on  the  right  side  of  the  second  equation  (3 .6. 1 .27)  becomes  CT Qr .  Much  can  be  said 

regarding  whether  these  equations  are  well  posed  and  stable.  From  Computational  Methods 
for  Linear  Control  Systems  [58],  a  series  of  tests  can  be  performed  to  check  those  conditions. 

Letting  N T N  =  Q  in  the  first  equation  of  (3.6.1.26)  or  Q  can  be  computed  from  N;  and 

letting  VTV  =  -  Qr ,  from  the  second  equation,  then  if  the  pair  (A,B)  is  controllable  and 

the  pairs  (N,  A)  and  (V,  A)  are  observable,  (i.e.,  O  =  [f  VA  VA2  •  •  •  VAn~x  ]  has 

the  full  rank  n  and  OtO  is  non-singular ),  the  Riccati  equations  must  have  unique  positive 
definite  solutions.  Thus,  there  exists  a  unique  optimal  control  u  *  (0)  which  minimizes 

equation  (3.6. 1 .3)  and  may  be  expressed  as  a  linear  state  feedback.  It  is  understood  that  A 
must  have  constant  coefficients  in  order  to  apply  the  LQR  theory. 

Boundary  conditions  imply  the  following: 

p*(0f)=  Sx*(Of)-  Sr(0f)=  k(0f)x*(0f)+  m(0f)  (3.6.1.28) 

V  x*(0f)  and  K^/)  ^en  k(0f)=  S  and  m(0f)=  - CTSr(0f )  (3.6.1.29) 

m(0  )  can  be  determined  by  integrating  backward  the  closed-loop  adjoint  system,  equation 
(3.6.1.29).  Then  m( 0)  is  known.  During  the  actual  control  run,  m(0)  is  used  in  the  forward 
equation,  the  second  equation  of  (3.6. 1 .27).  This  method  however,  will  not  be  used,  but  the 
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Riccati  equation  will  be  allowed  to  reach  a  steady-state  solution.  Therefore,  the  final 

boundary  conditions  will  not  be  of  consideration. 

3.6.2  Solution  of  the  LQR  based  on  the  Tschauner-Hempel  Equations  of 
Motion 

Proceeding  to  place  the  Tschauner-Hempel  equations  into  state  space  form,  these 
equations  can  be  represented  as: 


r-2 

ij"  +  2f'  =  an  0-6.2A) 

("  +  t=a( 


where  T],and  £  are  non-dimensionalized  coordinates  centered  at  the  target  or 

reference  spacecraft.  £  describes  cross  track  (outward  radial)  motion,  7  -  along  track,  and  £ 

out  of  the  nominal  orbit  plane  of  the  target  spacecraft.  Equation  (3.6.2.1)  can  then  be 
converted  into  state  variable  format,  as: 


r  =  2  7'+^+0{ 
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Upon  interchanging  rows  so  the  first  three  variables  are  position  and  the  second  three  are 
velocity,  then  equation  (3.6.30)  in  common  state  space  representation,  is: 
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3.6.2.2(b) 


The  nonlinear  term  in  this  matrix  can  be  adjusted  in  a  number  of  ways: 

1 .  When  it  is  assumed  that  r  remains  constant  (i.e.  equal  to  r(0)  =  —  ),  true  for  a  circle 


and  relatively  short  displacements,  then  the  term  becomes  equal  to  3. 

2.  If  simulation  is  started  at  perigee  or  apogee  then  evaluate  r  at  perigee  or  apogee 
respectively,  and  treat  as  constant  for  sufficiently  short  time  thereafter. 

3 .  If  several  orbits  are  needed  to  correct  the  disturbance,  then  use  an  average  value  for  r 
with 


h  = rx  v  = 


h  is  the  angular  momentum,  a  the  semimajor  axis  and  b 


is  the  semiminor  axis. 

4.  A  final  consideration,  is  to  update  1  and  2  in  a  piecewise  adapted  manner  along  the 
orbit. 

Although  the  T-H  equations  simplify  the  fundamental  matrix,  the  right  hand  side 
(RHS)  can  become  quite  cumbersome,  depending  on  the  type/location  of  actuators.  The 
terms  on  the  RHS  of  (3.6.2. 1)  contain  all  the  relative  accelerations  (i.e.  drag,  thrust  and 
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Earth’s  oblateness).  However,  since  the  mass  of  each  spacecraft  is  equal,  and  if  the  ballistic 
coefficients  are  the  same,  then  the  relative  drag  is  approximately  equal  to  zero.  Therefore 
its  contribution  will  not  be  considered  as  part  of  the  acceleration  term,  but  only  the  thrust  and 
oblateness  terms.  If  the  variables  cl,  c2,  and  c3  are  selected  to  include  thrust  and 

oblateness,  the  control  will  appear  as  uc  =  u-  d .  The  oblateness  term  will  be  included 
separately  and  the  development  will  proceed  as  below. 


The  contribution  of  oblateness  is  as  follows: 

G 


a  =  J\~y)  [Pk^OS<p)ir  -  p'k(C0S<p)i2 ]  (3. 6.2.3) 


tgq  is  the  equatorial  radius.  In  X,y,Z  coordinatesfto  be  normalized  in  the  r|,  C,  notation), 

where  u  =  co  +  0,  and  from  Battin  [18,52],  if  only  the  second  harmonic  term  is  considered , 
then  the  acceleration  terms  are  as  follows: 
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"0  0  0> 

0  0  0 

0 

0  0  0 

and  d  - 

0 

1  0  0 

aj 

0  1  0 

0 

^0  0  lj 

^  0y 

(3.62.5) 


where 


aj  = 


-3J1req1{\  +  ecosO)A 
a\\-e2)4 


Now  enough  information  is  known  to  execute  the  program  for  LQR.  This  LQR  program  will 
be  based  on  the  linearize  Tschauner  Hempel  equations  of  motion  with  respect  to  the  LVLH 
coordinate  frame.  The  origin  of  the  LVLH  is  the  mother  spacecraft  point  on  the  reference 
ellipse.  Again  the  t,  axis  of  the  frame  points  radially  outward  along  the  local  vertical 
direction,  the  r\  axis  is  along  the  direction  of  motion  on  the  reference  orbit,  while  the  C,  axis 

is  normal  to  the  reference  orbit  plane.  The  initial  state  vector  is  tj, £,£’,?]'  ,£')  and  the 

initial  perturbation  variable  as  x0  =  (8^,8rjy8^,^' ,f]' ,1^')  where 

8^-^  -  %ref  8jj  =  tj  -  J]ref  and  5£  -  £  -  £ref  ,  etc.  The  Tschauner-Hempel 

equations  will  be  linearized  only  once,  in  that,  station  keeping  will  be  considered  for 
spacecraft  along  track,  thus  the  reference  spacecraft,  can  be  taken  as  the  unperturbed 
position  of  any  daughter’s  position.(  If  some  other  configuration  is  to  be  considered,  such  as 
a  triangle,  then  these  equations  would  need  to  be  re-linearized  around  various  places  in  the 
orbit,  to  accommodate  an  inherent  two  dimensional  configuration).  For  example,  a  501  km 
(mother  -  daughter)  separation  from  a  nominal  500  km,  implies  a  one  kilometer  as  input,  as 
the  variation  of  the  daughter  spacecraft  from  its  unperturbed  position  along  track. 


67 


The  aj  term  in  equation  (3. 6.2.5),  and  others  when  necessary,  must  be  put  into  non- 

dimensional  form;  this  will  be  accomplished  in  the  coding,  being  mindful  of  theta 
dependance. 

Once  the  state  matrices  are  specified,  then  enough  information  is  available  to  solve 
the  differential  Riccati  equations.  A  couple  of  methods  exists  to  solve  equations  (3.6.1.26 
and  3.6.1 .27).  Since  in  general  k  and  m  are  an  n  x  n,  symmetric  matrices  (n  =  6 ),  implies 
n(n  +  l)/2  first-  order  differential  equations  must  be  solved.  These  equations  can  be 
integrated  numerically,  starting  at  0  =  0f  and  proceeding  backward  to  0  =  0O ;  k(0)  and 

m(0)  are  stored,  and  the  feedback  gain  matrix  is  determined  from  equation  (3.6.1.23). 

Now  most  Runge-Kutta  integration  routines  run  forward  in  time.  To  accomplish  the 
backward  integration,  the  RHS  of  (3.6.1.26)  is  multiplied  by  a  ( -1  ),  and  then  integrated 
forward  from  0  =  0O  (corresponding  to  t  =  0 ).  The  resulting  solution  is  reversed  and  shifted 
to  0  =  0f  (corresponding  to  t  =  T ).  The  second  part  of  the  simulation  is  to  find  the  optimal 
control  law  and  update  the  system  dynamics  by  integrating  forward  in  time.  Alternatively, 
to  determine  the  k  and  m  matrices,  for  an  infinite  time  process,  either  perform  a 
backward  recursion  integration  until  a  steady  state  solution  is  obtained  or  simply  solve  the 
nonlinear  algebraic  Riccati  equation  (ARE),  obtained  by  setting  k'  and  aw'  equal  to  0.  The 
second  part  of  the  simulation  still  is  to  find  the  optimal  control  law  and  update  the  system 
dynamics  by  integrating  forward  in  time.  This  method  will  be  chosen  to  solve  the  Riccati 
equation,  of  which  k  and  m  are  solutions.  The  program  contains  the  following  dimensions 

for  the  various  matrices:  A  e  R* 5x6  B  e  R{ 6x3  S  e  R6*6  R  €  jR3x3  Q  e  R6*6 
;  WiihRN  e  R6xl  and  x  e  R6xl  would  imply  a  calculated  value  of  control  u  e  RM , 

feedback  and  feedforward  respectively,  k  e  Rix6  &  me  RM  .  If  infinite  time  only  is 

considered  then  the  weighting  matrix  S  =  0,(  as  in  Eq.  (3.6. 1.3)),  since  only  an  engineering 
approximation  is  assumed.  The  MATLAB  function  [K  S,  e]  =  lqr(A,  B,  Q,  R,  N)  will  be 
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used  to  solve  the  ARE,  for  a  particular  Q  and  R;  N  is  defaulted  to  zero .  Further  descriptions 
show  that  this  function  also  calculates  the  optimal  gain  matrix  K  such  that  the  state-feedback 
law  u  -  -  Kx  (3.62.6) 

minimizes  the  quadratic  cost  function 

co 

J(u)=  \{xrQx+  uT Ru+  2xt  Nujdt  (3.62.7) 

o 

for  the  continuous  time  state  space  model  ~  +  d  (3  5  2. 8) 

The  oblateness  is  included  as  a  disturbance  term  in  (3.62.8).  The  state  feedback  gain  K,  Iqr 
returns  the  solution  S  of  the  associated  Riccati  equation.  As  a  note,  be  mindful  of  change  of 
variables.  For  instance,  in  the  MATLAB  routine  S  is  not  the  same  as  the  weighting  function 
S,  and  is  the  solution  to  the  ARE. 

^r5+  SA-(SB  +  N)R~l(BTS  +  Nt)+Q=  0  (3.62.9) 

and  the  closed  -  loop  eigenvalues  e  =  eig(A  -  B  *  AT).  (3.6.2.10) 

Note  that  K  is  derived  from  S  by  K  =  R~l(BTS  +  NT)  (3.62.1 1) 

or  u  =  -R~\BtS+  NT)x  (3.6.2.12) 

Again  the  problem  data  must  satisfy  the  following  limitations: 

•  (A,B)  stabilizable 

•  R>  0  and  Q-NR~'Nt>0  (3.6.2.13) 

•  {Q-  NR~XN\  A-  BR~lNT) 

The  last  expression  has  no  unobservable  pole  on  the  imaginary  axis.  Nevertheless, 
observability  is  not  being  considered  here.  That  is,  till  states  are  assumed  to  be  immediately 
available.  Furthermore,  for  the  purpose  of  this  research,  this  problem  is  assumed  to  be  in  the 


absence  of  noise.  ARE,  equations  (3.6.1.27),  will  be  solved  using  MATLAB  notation  and 
variables.  The  first  of  these  is  the  same  as  equation  (3.6.2.9)  with  N  =  0.  The  second  has  the 

solution  m=(AT-  SBR^BTy\CTQr -  Sd)  (3.6.2.14) 

When  the  gain  reaches  a  steady-state  value  and  the  closed-loop  plant  is  stable  then  the 
optimal  tracker  is  given  by 


-  m'  =  (A  -  BK( oo  ))Tm  +  CTQr  +  Sd 

u=-K(n)x+  R'lBT(AT-SBR-'BTylCTQr- R-'BT(Ar  -  SBR~}BTy'Sd  (3-6.2.15) 


The  second  equation  of  equation  (3.6.2. 1 5),  represents  the  control  needed.  The  first  term  of 
this  same  equation,  represents  the  feedback,  second  feedforward  or  tracking  and  the  third 
allows  for  disturbance.  The  form  of  the  control  is  similar  to  that  which  would  result  had  the 
disturbance  been  included  as  part  of  the  B  matrix. 

Letting  the  steady  state  track  and  disturbance  be  represented  as, 

track  =  R-'Bt(At  -  SBR~'BTy'CTQr 

and  (3.6.2.16) 

disturbl  =  R~'BT(AT  -  SBRT'Bry'Sd 


the  state  becomes 


x'  =  Ax+  Bu+  d  -  Ax+  B{-  K(<x>  )x  +  track  -  disturbl)  +  d 


(3.6.2.17) 
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Upon  substitution  for  u,  this  equation  can  be  written  as 


x'  =  (A  -  BK)x+  Btrack-  Bdisturb\+  d  =  Ax+  fJ  (3.6.2.18) 
which  has  the  solution  of  the  form 

°i 

x  =  exp(A  (6f  -  0o))jco  +  exp(A  6f  )  J exp(A  0 )  fJdO  (3.6.2.19) 

*0 


where  A  =  (A  -  BK) 

The  auxiliary  signal  (the  second  part  of  equation  (3.6.2.19))  can  be  integrated  to  be 


/=  J(expA(tf/-^))/y^=(/'/AXl-expA(^-«l,))  (3.6.2.20) 


In  the  infinite-horizon  limit,  this  signal  is  a  constant  and  is  represented  by 

/( oo )  =  /  A  =  A  _1  ( Btrack  -  BdisturbX  +  d)  (3.6.2.21) 

The  solutions  to  equations  (3.6.2.15)  and  (3.6.2.19)  are  easily  coded  in  the  program  strjqr. 
This  code  allows  one  the  option:  (l).the  use  of  feedforward  (tracking)  or  (2)  not  to  use 
feedforward  (tracking). 

Function  lqr  can  be  used  alone  or  in  conjuction  with  a  general  calculation  of  the 
optimal  feedback  for  constant  gain,  first  presented  by  Moerder  and  Calise  [59]  based  on  the 
Lyapunov  equation.  The  Lyapunov  equation  is  a  variation  of  the  Riccatti  equation  obtained 
by  making  a  transformation  of  the  fundamental  matrix.  The  Lyapunov  equation  may  be 
solved  using  the  subroutine  lyap.m  in  MATLAB  (Control  System  Toolbox).The  initialization 
of  the  program  requires  an  initial  stabilizing  output  feedback  kO,  ( this  kO  can  be  provided 
by  lqr),  the  maximum  number  of  iterations  N,  the  magnification  a  of  Ak,  and  the  tolerance. 
The  weighting  Q  and  R  are  also  inputs  to  the  program.  The  number  of  iterations,  the 
tolerance,  and  the  magnification  can  determine  the  successful  termination  of  the  program. 
If  a  is  “large”  the  matrix  ‘A  -  BKC’  may  not  be  Hurwitz.  That  is,  it  may  not  contain 
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negative  real  parts  of  the  eigenvalues.  In  such  cases,  the  algorithm  will  terminate  with  a 
warning.  To  alleviate  the  problem,  decrease  a.  If  the  tolerance  is  very  “small”  then  the 
program  may  terminate  before  achieving  the  desired  tolerance.  In  that  case,  either  increase 
N  or  the  tolerance. 

In  the  interest  of  time,  that  is  for  a  timely  completion  of  this  phase  of  the  research, 
lqr  will  be  used  in  isolation  of  lyap.m.  Thus  any  reference  to  this  function  will  be 
‘commented  out’.  Refining  of  the  optimal  software  can  be  explored  in  future  research.  The 
advantage  of  the  combination  software  usage  is  that  quicker  satisfying  weighting  matrices 
may  be  determined. 

As  a  final  comment  regarding  the  criterion  functional,  the  system  is  desired  to  be 
optimal  but  must  be  very  exact  in  the  sense  in  which  the  system  is  optimal.  Mathematical 
expressions  must  be  found  to  measure  how  the  system  must  be  optimal  in  comparison  with 
other  system.  The  system  response  is  usually  of  most  interest.  The  exact  form  of  S  and  Q 
is  to  be  fixed  by  the  designer  at  the  outset.  Thus  a  number  is  assigned  to  the  response 
obtained  by  each  control  law  u,  and  the  optimum  system  is  that  whose  control  law  gives  the 
minimum  cost  function.  The  choice  of  Q  is  dictated  by  the  relative  importance  of  each  state 
over  some  time  interval  or  [0O  <  0  <  0f  ]  in  this  case.  Fast  decay  implies  large  control 

or  large  Q  implies  large  u.  The  choice  of  S  is  dictated  by  the  relative  importance  of  each 
state  at  the  final  time  or  designation.  The  relative  magnitudes  of  Q  and  R  are  in  proportion 
to  the  relative  values  of  the  response  and  control  energy.  The  larger  Q  is  relative  to  R,  the 
quicker  the  response  and  the  higher  the  gain  of  the  system. 

The  calculated  thrust  is  a  non-dimensional  theta  dependent  quantity.  Carter  and 
Brient  [29],  stated  that  the  Tschauner-  Hempel  vector  x(0)  described  by  the  Tschauner- 

Hempel  equations  can  be  obtained  from  the  transformation  x(0)  =  r(0)X(0)  (3.6.46) 

where  X(0)  is  the  actual  relative  position  of  the  spacecraft  in  terms  of  0.  If  the 
solution  vector  to  these  equations  is  used,  then  this  would  imply  that 
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(3.6.2.22) 


v(0)  =  r(0)X'(0)  +  r'(0)X(0) 

m-rwxv) 

=>  X  {0}  -  ... 


v(0  = 


dX\0)d0 

dtdO 


X'(0)0 


Similarly,  and  using  the  fact  that 


it  can  be  shown  that 

4 

0(0=  ^-{2rJWr'W[vOT-  r'(e)X(0)]} 

4 

+  r"(9)JT(tf)]-  2r' (O)r1  (0)(v(O) - 


using  again  the  approximation  r(0)  = 


h 2 


therefore 

a(t)  =  ^-{2 /(#)-r'(W)f 
+/»{r(0)[«(0)-  c”(W)]- 2r'(0)(v(0) 


From  equation  (3.6.2.22) 


(3.6.2.23), 


(3.6.2.24) 

r\e)X(0))} 


(3.6.2.25) 


(3.6.2.26) 

r’{9)X{0j)} 
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a{0 )  =  2r'X'  +  rX"  +  *  2r'JT  +  r"X 


(3.6.2.27) 


then 


a(/)  =  ^-{y(«)[vW-r'WJfW]} 

+n{r(0)[2r'(9)X'(0)\-  2r'(0)(v(0)-  r’(9)X(0))} 


If  equation  (3.6.2.25)  is  not  used,  but  instead 


(3.6.2.29) 


then 


a(t)  = 


{ 2r2r'(v-r’X)+2r4rX ' 


r'r3(v-  rX)\ 


(3.6.2.30) 


Equation  (3.6.2.28)  or  (3.6.2.30 )  are  also  programed  using  MATLAB  and  is  included  in  the 
program  str  lqr. 

Assuming  all  states  are  immediately  available  and  in  the  absence  of  noise,  a 
parametric  study  will  be  investigated  using  various  weighting  functions.  The  task  then  is  to 

select  the  combination  of  values  for  the  state  and  control  weighting  matrices  which  will  best 

answer  questions  for  a  given  scenario,  while  maintaining  other  mathematical  behaviors  such 
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as  positive  definiteness  and  positive  semi-definiteness(  for  purposes  of  stability).  The  three 
relevant  questions  the  LQR  should  address  are  (1)  How  much  control?  (2)  How  much 
overshoot?  and  (3)  How  much  time  is  needed  for  convergence?  All  of  which  equate  to 
values  needed  to  remove  the  error.  Although  the  designed  matrices  were  given  an  out-of- 
plane  ‘z’  component  as  input,  it  was  not  until  Figures  3.6.2.6  and  3.6.2.7  that  a  ‘z’ 
component  error  was  inputted.  Thus  there  are  no  ‘z’  responses  for  the  earlier  Figures  since 
a  z  response  at  this  point  is  irrelevant  for  a  coplanar  model.  A  discussion  per  Figure  is 
delineated  separately.  The  preliminary  results  shown  here  assume  the  initial  LQR  correction 
begins  near  perigee,  at  a  true  anomaly  angle  of  45  °.  The  periodic  system  matrix  appearing 
in  the  Tschauner-Hempel  equations  is  evaluated  at  that  true  anomaly  angle.  If  the  responses 
occur  in  a  relatively  short  time,  it  is  assumed  that  this  value  could  be  used  throughout  the 
maneuver.  For  longer  time  responses  this  matrix  would  have  to  be  re-evaluated  in  a  piece- 
wise  adaptive  manner.  For  all  results  presented  here,  the  quadratic  cost  functional  with  the 
infinite  time  upper  bound  and  N  =  0  as  in  Eq.  (3.6.2.7)  has  been  selected  as  the  quadratic  cost 
function. 
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R  Dag  =1  1  1 


Q=1  1  1  1  1  1 


time  (sec)D  time  (sec)D 


Figure  3.6.2. 1 

Figure  3.6.2. 1  will  be  considered  as  the  reference  Figure  for  the  others.  The  diagonal 
elements  of  the  control  R  and  state  Q  weighting  matrices  are  given  the  value  of  ones.  The 
off  diagonal  terms  are  assumed  to  be  zero.  The  state  here  is  defined  as 
(x  y  z  x  y  z)  where  x  is  in  the  radial  cross-track  direction  in-plane,  y  is  in  the 

along  track  direction,  and  z  represents  out-of-plane  displacement.  The  center  of  this  moving 
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coordinate  system  is  taken  at  the  nominal  position  of  any  daughter  satellite  in  the  orbit.  Thus 
a  displacement  of  1  km  along  track  actually  represents  a  displacement  of  dmotherKUughttr  + 1  km 
along  track,  where  d  mothCT.<jaUghter  is  the  desired  mother-daughter  separation  distance  for  this 
application. 


RQag  =3  3  3 


Q  =1  11111 
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Figure  3.6.2.2 


The  diagonal  elements  of  R  are  increased  to  3.  As  such,  the  control  effort  decreases  as 
compared  with  Figure  3 .6.2. 1 .  Also  some  of  the  responses  (velocity  and  acceleration)  appear 
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to  be  slightly  more  sluggish  or  slow  to  respond. 


RQag  =1  1  1 
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Figure  3. 6.2.3 


The  penalty  on  all  state  components  is  increased.  Some  of  the  overshoots  are  noticeably 
reduced  as  compared  with  Figure  3.6.2. 1  responses.  The  responses  of  the  x  component  of 
acceleration  and  position  are  much  faster.  The  peak  control  effort  components  are  an  order 
of  magnitude  larger  than  in  Figure  3.6.2. 1. 


78 


R  Dag  =111 


Q=1  50  1  1  50  1 


time  (sec)o  time(sec)D 


Figure  3. 6.2.4 


Now  penalties  on  the  y  component  of  the  state  matrix  and  also  on  the  y  component  of  the 


state  matrix  are  increased  50  fold.  The  end  result  is  almost  an  order  of  magnitude  increase 
in  control  effort,  and  a  slight  improvement  in  reducing  some  of  the  transient  overshoots  as 
compared  with  the  results  in  Figure  3.6.2. 1. 
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Figure  3.6.2.5 

In  this  example,  the  split  weighting  strategy  is  employed.  R  is  penalized  more  than  before, 
requiring  in  some  cases  almost  4000  sec  (twice  the  time  of  Figure  3.6.2.1)  to  reach  steady 
state;  this  results  in  an  order  of  magnitude  reduction  in  control  effort  u. .  Transient  responses 
are  noticeably  more  sluggish. 
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Figure  3. 6.2.6 

Although  the  previous  results  emphasized  coplanar  type  maneuvers,  here  an  out-of-plane 
error  in  position  of  1  km  is  now  given  at  the  45 0  true  anomaly  point,  in  addition  to  the  1  km 
error  along  track.  The  weighting  functions  are  the  same  as  those  of  Figure  3.6.2. 1.  This 
Figure  depicts  consistence  as  it  relates  to  input  error.  Furthermore,  damped  simple  harmonic 
motion  (SHM)  is  demonstrated  as  predicted  by  the  out-of-plane  Tschauner-Hempel  equation. 
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km 

Figure  3.62.1 

Figure  3.6.2.7  is  a  three  dimensional  depiction  of  Figure  3.6.2.6.  The  weighting  elements 
selected  here  by  no  means  infer  that  these  are  the  optimal  choices  or  combination  of  choices, 
only  to  demonstrate  capability  and  workability. 

Preliminary  results  assume  the  initial  LQR  correction  begins  near  perigee,  at  a  true 
anomaly  angle  of  45  °.  The  preliminary  results  here  validate  the  controllability  and  stability 
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of  the  LQR  control  strategy  and  verify  convergence  in  a  relatively  short  time  interval  in 
response  to  an  initial  displacement.  The  response  time  is  much  shorter  than  that  of  Tan,  et. 
al.[38]  or  Schaub,  et.  al.[51].  Both  used  a  Lyapunov  type  of  approach.  Tan  applied  the 
approach  to  the  osculating  orbital  elements,  while  Schaub  used  the  mean  orbital  elements. 

3.6 3  Modification  to  Station  Keeping  Strategy  due  to  the  Earth’s  Planetary  and 
Lunar  Perturbations 

To  understand  the  effects  of  Earth’s  planetary,  and  Lunar  perturbations,  let’s  begin 
with  a  description  of  the  software  and  its  inputs.  The  Ada  Simulation  Development  System 
(ASDS)[21]  or,  more  commonly,  BG14,  has  the  following  environment  models  available: 
Atmosphere,  Gravity,  Gravity  Gradient  Torque,  Third  Body  Effects,  Tidal  Forces  due  to  Sun 
and  Moon,  Solar  Radiation  Pressure  and  Solar  System  Models,  including  ephemeris  and 
planet  physical  parameters.  Two  ephemeris  modeling  capabilities  exist  in  ASDS.  One  uses 
an  analytical  approach  to  computing  planetary  ephemerides  (V an  Flandem),  while  the  other 
one  uses  actual  Jet  Propulsion  Laboratory  (JPL)  data  for  the  planets  and  satellites.  The  Van 
Flandem  model  gives  low-precision  (1  degree)  formulas  for  geocentric  and  heliocentric 
positions  of  the  Sun,  Moon,  and  planets,  which  are  valid  for  any  epoch  within  300  years  of 
the  present.  In  the  JPL  version,  an  ephemeris  of  the  Moon  and  nine  planets  has  been 
numerically  integrated  from  1411  BC  to  3002  AD.  A  long  ephemeris  has  utility  for 
comparison  with  both  historical  observations  and  analytical  theories.  The  environment  data 
set  provides  inputs  as  central  body  and  the  perturbing  forces.  It  is  seen  that  Central_  Body 
=  Earth,  Third_Body_Plantes=  Sim  and  Moon,  Gravity _Model=wgs_84  (forces  the  use  of 
an  Earth  fixed  system),  Gravity(with  various  order  and/or  degree)  and  Solar_Radiation  is  set 
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to  either  ‘on’  or  ‘off.  Note  that  the  solar  system  model(ephemeris)  used  is  the  JPL 
Ephemeris  unless  the  date  specified  lies  outside  the  valid  ephemeris  dates,  in  which  case  Van 
Flandem  is  used. 

Higher  harmonics  of  perturbations  in  three  dimensions,  can  be  seen  from  the  more 
general  form  for  the  Earth’s  gravitational  potential  which  can  be  written  as 
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Jk  is  defined  as  Jk  =  -C°k  .Values  of  the  various  coefficients  are  obtained  from  satellite 

observations.  The  k’s  correspond  to  degrees,  whereby  the  j’s  depict  the  order.  As  the  degree 
and  order  increase,  the  order  of  magnitude  effect  of  each  harmonic  on  satellite  motion 
decreases. 
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Figure  3.6.3.1(a)  30  Day  Time  History  of  Separation  Distance 

between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 

From  Figure  3.6.3.1(a),  the  degree  as  well  as  order  was  set  to  equal  20.  These  higher 
harmonics  seemed  to  have  negligible  effect  in  comparison  with  Fig  3.4.1(b),  for  a  satellite 
with  an  initial  separation  of  400  km  for  a  30  day  orbit. 

On  the  other  hand,  higher  harmonics  do  contribute,  although  slightly,  for  orbits  in 
excess  of  1 00  days.  Here  the  question  of  collision  is  answered;  that  is,  if  no  other  corrections 
are  given,  “what  happens?”  The  effects  of  including  these  higher  harmonics  (20  by  20) 
delays  the  collision  time  by  25  days,  otherwise  collision  time  would  have  been  reached  at 
approximately  125  days,  instead  of  the  now  150  days.  See  Figures  3.6.3. 1  (b  and  c). 
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Figure  3.6.3.1(b)  1000  Day  Time  History  of  Separation  Distance 

between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 
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Figure  3.6.3.1(c)  1 000  Day  Time  History  of  Separation  Distance 

between  Mother  and  Daughter  Satellite  Non-Keplerian  Orbit 

In  summary,  for  relative  short  orbit  periods  (i.e.  30  days  =  236+  orbits),  the  Earth’s 
planetary  and  Lunar  perturbation  exhibited  negligible  effect.  For  higher  harmonics,  the 
Earth’s  planetary  and  Lunar  perturbation  became  noticeable  for  orbits  of  100  days  or  more. 
It  is  also  shown  that  higher  orders  and  degrees  of  perturbation  increased  the  lifetime  of  the 
satellite  by  prolonging  the  time  before  zero  separation  is  reached.  The  shifting  of  lines  of 
apsides  for  this  elliptical  orbit,  showed  an  order  of  magnitude  improvement  in  reaching  the 
time  of  collision,  as  compared  to  that  of  Badesha’s  circular  orbit  [32]. 
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IV  THE  IMPLEMENTATION  OF  MAINTAINING  CONSTANT 


DISTANCE  BETWEEN  SATELLITES  IN  ELLIPTIC  ORBITS 

4.1  INTRODUCTION 

The  scientific  objectives  of  the  Earth  observation  program  are  becoming  more 
autonomous  and  more  ambitious;  this  has  created  needs  for  innovative  technical  approaches 
to  achieving  and  maintaining  constellation  and  formation  flights  of  spacecraft.  The  trend  to 
develop  small  low-cost  spacecraft  has  led  many  to  recognize  the  advantage  of  flying  multiple 
spacecraft  in  formation  to  achieve  the  correlated  instrumentation  formerly  possible  only  by 
placing  many  instruments  on  a  single  platform. 

A  study  was  conducted  of  proposed  NASA  and  ESA  constellation  configurations 
which  would  measure  and  study  upper  atmospheric  phenomena.  The  Auroral  Cluster 
(Multiscale)  System,  the  Distance  Measurement  System,  the  Orbiting  Interferometer  System, 
as  suggested  by  NASA  for  LEO  missions  together  with  the  Solar  Stereo  System  in 
heliocentric  orbit  were  all  considered  as  possible  baseline  or  “strawman”  configurations.  In 
addition  information  was  also  obtained  from  the  ESA  web  page  on  the  proposed  ESA  Cluster 
mission  with  the  objective  of  determining  physical  processes  involved  in  the  interaction 
between  the  solar  wind  and  the  magnetosphere  by  visiting  key  regions  such  as  the  polar  cusps 
and  magnetotail. 

After  reviewing  the  candidate  configurations  it  was  decided  to  select  the  Auroral 
Multiscale  Mission  (AMM)  as  a  baseline  or  “strawman”  model  for  this  research:  a  brief 
description  and  sketch  of  this  system  is  given  in  Fig.4. 1  and  Table  4. 1 .  Initially  for  this  study 
the  “strawman”  configuration  would  be  based  on  four  spacecraft  in  the  same  plane 
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SingleAMM  Spacecraft  /  Stowed  Configuration 


Launch  Configuration 


Fig.  4.1  Auroral  Multiscale  Mission  System 
Table  4.1.  Strawman  Configuration 


Mission 

Multi-point  3 -dimensional  data  collection  of  auroral  phenomena 

Orbit 

600x8000  km,  83°  inclination 

Launch  Vehicle 

Taurus  (2110  vehicle),  Insertion  Mass  of  465.5  kg,  argp=203.1° 

Spacecraft  Size 

40  inch  diameter  (across  flats) 

Spacecraft  Mass 

90  kg  each 

Science  Payload 

Ion/Electron  spectrometer,  UV  Auroral  Imager,  Magnetometer,  Electric 

Fields  (3-axis) 

Position  Knowledge 

GPS,  Knowledge  to  1 00  meters  (3  sigma) 

Attitude  Knowledge 

0.01°  Star  Tracker  (referenced  to  magnetometer),  star-tracker  implementa¬ 
tion.  Spinning  sun-sensor  /  magnetometer  provides  coarse  attitude. 

Power 

Solar  array  capability:  40  watts  orbit  average  power 

Energy  Storage:  Dual  IP  ACS  Flywheel  momentum  bias  system  used  for  both 
power  and  attitude  control 

principally  along  the  orbit  track. 
The  perigee  altitude  is  selected  to 
be  600  km,  the  apogee  altitude  is 
selected  to  be  8000  km.  and  the 
nominal  separation  distance 
between  two  adjacent  satellites  is 
taken  to  be  500  km  Without  any 
perturbation  or  any  control,  the 


Fig. 4. 2  Variable  Velocities  &  Phase  Distances 

velocity  at  perigee  and  apogee  as 
well  as  the  separation  distances  are  shown  in  Fig.4.2.  From  Fig.4.2  it  is  seen  that  the 
separation  distances  at  perigee  are  more  than  twice  the  separation  distances  at  apogee.  To 
maintain  constant  separation  distance  in  such  an  orbit,  it  would  be  necessary  to  correct  the 
orbit  continuously;  the  tremendous  amount  of  energy  required  makes  this  strategy  unfeasible. 

A  novel  idea  for  implementation  of  constant  separation  distances  for  the  AMM 
mission  is  developed49.  The  four  satellites  are  launched  by  a  single  vehicle;  the  method  to 
distribute  the  satellites  into  their 
positions  using  the  least  maneuver 
energy  is  given  in  this  chapter. 

Because  of  various  perturbations 
(mainly  caused  by  the  J2  effect)  and/or 
initial  errors,  the  separation  can  not  be 
maintained  without  control.  There  are 
several  mathematical  models  that  can  be  Fig. 4. 3  Angle  Between  Orbits  of  Two  Satellites 
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used  to  design  the  controllers.  The  Clohessy- Wiltshire  equations22  (sometimes  referred  to  as 
the  Hill-Clohessy- Wiltshire  equations60 )  assume  circular  orbits  and  a  spherical  Earth  model, 
and  the  equations  of  motion  of  the  orbiting  spacecraft  are  linearized  relative  to  the  rotating 
frame  of  the  reference  spacecraft.  Here,  however,  the  orbit  is  elliptic. 

In  this  chapter  by  controlling  the  osculating  orbital  elements  the  constant  separation 
is  successfully  maintained.  The  advantage  is  taken  of  celestial  mechanics  insight  to  avoid  the 
correction  of  orbital  elements  at  ill-suited  times.  Since  the  model  involving  the  osculating 
orbital  elements  is  inherently  nonlinear  and  time-variant,  control  logic  based  on  a  Lyapunov 
function  is  applied.  The  control  is  much  simpler  than  the  one  of  Schaub  et  al51,  therefore  it 
is  easier  to  implement  in  engineering  practice.  All  simulations  are  performed  by  MATLAB 
and  the  BG14  orbital  propagator 21 . 

4.2  REVIEW  OF  STRATEGY  FOR  MAINTAINING  DISTANCE  IN  ELLIPTIC  ORBITS 

In  order  to  maintain  a  constant  distance  between  two  satellites  in  elliptic  orbits,  the 
semi-major  axis  of  the  orbit  of  the  second  satellite  should  be  shifted  by  a  very  small  angle 
(  something  like  1.37° )  with  respect  to  the  orbit  of  the  mother  satellite  in  order  to  achieve 
the  constant  separation  distance.  Both  orbits  are  in  the  same  plane  (Fig.4.3).  If  there  are  two 
satellites  flying  in  the  orbits  shown  in  Fig.4.3,  both  in  the  counter  clockwise  direction,  at 
exactly  the  same  time  that  they  arrive  at  their  perigees  and  apogees,  respectively,  the  distance 
between  the  two  satellites  is  defined  as  the  geometric  distance,  roughly  the  distance  caused 
by  the  shift  of  the  semi-major  axis  of  the  orbit;  if  the  satellites  are  flying  in  exactly  the  same 
orbit  (shown  in  Fig.4.2),  the  distance  between  the  adjacent  satellites  is  defined  as  the  phase 
distance,  i.e .  the  distance  accounted  for  during  the  time  in  which  the  first  satellite  will  reach 
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the  second  satellite’s  current  position.  From  Fig.4.3,  it  is  seen  that  the  geometric  distance 
at  the  perigee  is  smaller  than  the  one  at  the  apogee.  The  strategy  here  is  that  if  the  phase 
distance  is  compensated  by  the  geometric  distance,  then  the  resulting  separation  distance 
between  the  adjacent  satellites  is  maintained  to  be  essentially  constant49.  The  details  of  the 
calculation  are  presented  as  follows: 

From  Fig.4.3,  if  the  angle  between  the  two  semi-major  axes  is  a,  the  distances 
between  satellites  1  and  2  at  perigee  and  apogee  (the  geometric  distances)  are 


a  cc 

2 (Re+perigee  altitude)sia(-)  and  2 (Re+ apogee  altitude) sin(— ) 


respectively,  where  Re  is  the  radius  of  the  Earth.  Let  the  phase  distances  at  perigee  and 
apogee  be  Pp  and  Pa,  respectively  (Fig.4.2),  We  try  to  make  the  resulting  distance  at  apogee 
the  same  as  the  one  at  perigee,  e.g.  512  km.  i.e. 


2  {Re+perigee 
2  (Re+apogee 


altitude)  sin(— )  +Pp 
2 


altitude)  sin( — )  -Pa 
2 


512 

512 


(1) 

(2) 


From  Fig.4.2,  it  is  noticed  that 


Pp_  _  517.28 
Pa  251.13 


(3) 


The  solution  of  equations  (l)-(3)  is  a  =1.37°,  Pp=344.705  km,  Pa=167.295  km.  That  means, 
satellite  1  arrives  at  its  own  perigee  (or  apogee)  39.5  seconds  later  than  satellite  2  does. 

An  example  of  the  constellation  in  an  orbit  described  above,  with  nominal  separation 
distance  of  500  km.,  is  given  here.  The  simulation  is  done  by  MATLAB  and  BG14;  the 
simulation  results  are  identical,  shown  in  Fig.  4.4,  which  shows  the  simulation  results  for  16 


92 


Keplerian  orbits.  The  maximum 


distance,  512  km,  occurs  at 
perigee  and  apogee;  the  minimum 
distance,  488  km,  takes  place 
between  the  perigee  and  apogee. 
Without  perturbations  and 
subsequent  control  this  distance 

0  2  4  6  8  10  12  14  16 

would  be  maintained  forever;  the 
Fig. 4. 4.  Separation  Distance,  Keplerian  Orbit  drift  from  the  nominal  value,  500 

km,  is  about  ±2.4%. 

4.3  INITIAL  SEPARATION  DEPLOYMENT 

It  is  assumed  that  the  four  satellites  are  launched  by  a  single  vehicle,  see  Fig.4.1; 
therefore,  it  is  important  to  study  the  method  to  separate  the  satellites  using  the  least 
maneuver  energy.  After  comparison  of  several  maneuver  methods,  the  following  approach 
was  selected,  because  all  the  work  done  by  the  thrusters  is  used  to  augment  the  energy  which 
is  required  from  the  transient  orbit  to  the  final  orbit. 

There  are  two  tasks  for  initial  separation:  one  is  to  cause  a  shift  in  the  angle  (about 
1.37° )  between  the  semi-major-axes  of  the  orbits  of  the  adjacent  satellites;  the  other  is  to 
create  a  phase  difference,  i.e.  according  to  reference  1,  the  second  satellite  arrives  at  its 
perigee  39.5  seconds  earlier  than  the  first  one  does.  It  is  well  known  that  if  a  satellite  travels 
in  a  circular  orbit,  its  speed  is  the  same  at  any  point;  and  if  the  satellite  is  accelerated  to  a 
larger  speed  (V)  along  the  track  at  some  point,  say  A,  then  the  satellite  will  go  into  an 
elliptical  orbit,  whose  perigee  is  A,  and  the  apogee  depends  on  the  speed  V.  Therefore,  the 
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four  satellites  are  launched  by  a  single  vehicle  into  a  circular  orbit  600  km  from  the  Earth’s 
surface  with  a  required  ascending  node  and  inclination  angle.  From  this  orbit,  if  any  satellite 
is  accelerated  to  a  larger  speed  along  the  track,  its  perigee  must  be  600  km.;  and  it  is  easy  to 
cause  a  shift  in  the  angle  between  the  semi-major-axes  of  the  orbits  of  the  adjacent  satellites. 

®  The  four  satellites  are  traveling  in  the  circular  orbit  (orbit  5  in  Fig.4.5),  the  speed 
of  this  orbit  is  7.558  km/s.  When  they  move  to  the  position  of  the  required  argument  of 
perigee  of  the  mother  satellite,  the  mother  satellite  will  be  released  and  accelerated  along  the 
track  to  8.770  km/s,  the  perigee  speed  of  the  final  orbit.  Thus,  the  mother  satellite  will  travel 
in  an  orbit  whose  perigee  is  at  this  point,  600  km  from  the  Earth’s  surface,  and  the  apogee 
is  8000  km  from  the  Earth’s  surface  (the  orbit  1  in  Fig.4.5).  The  other  three  daughter 


Fig.4.5  Deployment  of  the  Four  Satellites  in  Constellation 

the  unit  for  both  coordinate  axes  is  meter,  s’s  inside  the  figure  mean  seconds,  the  position  10981  s  is  identical  with  the  point 
0.  but  at  point  0.  all  four  satellites  are  joined  together,  at  10981  s  they  are  separated  500  km.  apart. 
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satellites  are  still  moving  in  orbit  5  at  this  moment. 

(D  22.56  seconds  after  the  release  of  the  mother  satellite,  the  three  daughter  satellites 
have  swept  through  1 .37°  in  the  circular  orbit.  At  this  moment,  if  the  first  daughter  satellite 
is  released  and  accelerated  along  the  track  also  to  8.770  km/s,  the  first  daughter  satellite  will 
travel  in  an  orbit  whose  perigee  is  at  this  point,  600  km  from  the  Earth’s  surface,  and  the 
semi-major-axis  is  1.37°  apart  from  the  one  of  the  mother  satellite.  At  this  moment  , 
however,  the  mother  satellite  already  left  its  perigee  22.56  seconds  earlier;  that  means  the 
phase  of  the  first  daughter  satellite  is  22.56  seconds  later  than  the  one  of  mother  satellite.  But 
it  is  required  that  the  phase  of  the  first  daughter  satellite  be  39.5  seconds  earlier  than  the  one 
of  the  mother  satellite.  Therefore  we  let  the  first  daughter  satellite  travel  in  a  transient  orbit 
whose  period  is  22.56+39.5=62.06  seconds  less  than  the  period  of  the  mother  satellite’s 
orbit.  The  period  of  the  mother  satellite’s  orbit  is 

r=27tv/aVY=10981.31s 

where  a  is  semi -major-axis  of  the  mother  satellite’s  orbit;  y,  the  gravitational  constant, 
=3.986005><1014,  so  the  period  of  the  first  daughter  satellite’s  transient  orbit  T,,  its  semi¬ 
major  axis  a,  and  the  speed  v,  at  its  perigee  are 

j  I  i 

r  =r-62.06  =10919.25  s  a ,=(—y2)3  =10637.867  km 

2n 

cl=al -(perigee +fle)  =  10637.867 -600 -6378. 137  =3659. 73  km. 

el=cl/al  =0.344028552  b\=a\-c\,  5^9988.523  km.  ^ 

2  , -  b.h. 

pl=blial  vi~ . . . — - =8.762  km/s 

where  R,,  is  the  radius  of  the  Earth.  So  the  first  daughter  satellite  is  accelerated  to  8.762  km/s 
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instead  of  8.77  km/s.  Thus  it  will  travel  in  a  smaller  orbit  (orbit  2  in  Fig.4.5),  whose  semi- 
major-axis  is  10637.867  km.,  and  period  is  10919.25  seconds.  That  means  the  first  daughter 
satellite  returns  to  its  perigee  in  62.06  seconds  less  than  the  mother  satellite,  but  the  first 
daughter  satellite  is  released  22.56  seconds  later  than  the  mother  satellite,  so  its  phase  is 
39.50  seconds  ahead  of  the  mother  satellite,  i.e.  the  phase  distance.  Once  the  first  daughter 
satellite  returns  to  its  perigee,  it  will  be  accelerated  again  until  the  speed  reaches  8.770  km/s. 
So  it  will  finally  travel  in  an  orbit  with  the  same  orbital  elements  as  the  mother  satellite,  i.e. 
the  semi -major-axis,  eccentricity,  inclination,  ascending  node,  except  for  the  the  argument 
of  the  perigee  (the  difference  is  1 .37°),  and  the  mean  anomaly  (the  difference  is  n><39.5). 

©  Similarly,  22.56  seconds  after  the  release  of  the  first  daughter  satellite,  the 
remaining  two  satellites  have  swept  through  another  1.37°  in  the  circular  orbit.  At  this 
moment  the  second  daughter  satellite  is  released  by  using  the  same  method  as  in  Eqs.  (4), 
i.e.,  O  T2=Tj-62.06,  ©  Calculating  a2  from  T2,  ©  Calculating  from  ©  Calculating  e2 
and  b2  from  a2  and  c  2  then  p2  and  h2,  @  Eventually  v2  =8.754  km/s.  So  the  second  daughter 
satellite  is  accelerated  along  the  track  to  8.754  km/s.  In  this  way  the  second  daughter  satellite 
will  travel  in  an  orbit  (orbit  3  in  Fig.4.5)  which  is  1.37°  away  from  the  first  daughter  satellite 
and  39.50  seconds  in  phase  ahead  of  the  first  daughter  satellite.  When  the  second  daughter 
satellite  returns  to  its  perigee,  it  will  also  be  accelerated  to  8.770  km/s. 

®  Similarly,  22.56  seconds  after  the  release  of  the  second  daughter  satellite,  the  last 
daughter  satellite  is  released  and  accelerated  along  the  track  to  8.746  km/s  so  that  the  last 
daughter  satellite  will  travel  in  orbit  4  in  Fig.4.5.  It  will  also  be  accelerated  to  8.770  km/s 
when  it  returns  to  its  perigee. 

The  simulation  results  are  shown  in  Fig.4.5,  where  the  unit  for  both  coordinate  axes 
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Table  4.2.  Total  Av’s  Required  for  the  Four  Satellites  and  the  Moments  when 
They  Should  be  Accelerated,  Respectively. 


THE  FIRST  MANEUVER 

THE  SECOND  MANEUVER 

Time, sec 

Av  required,  m/s 

Time,  sec. 

Av  required,  m/s 

Mother  Sat. 

0 

8770.1-7557.9=1212.2 

10981.31 

0 

Daughter  Sat.  1 

22.56 

8762.0-7557.9=1204.1 

10941.81 

8770-8762=8.06 

Daughter  Sat.2 

45.12 

8753.9-7557.9=1196.0 

10902.31 

8770-8754=16.2 

Daughter  Sat.3 

67.68 

8745.6-7557.9=1187.8 

10862.81 

8770-8746=24.4 

is  meter,  s’s  inside  the 
figure  mean  seconds, 
i.e.  the  time  for  release 
of  the  mother  satellite 
is  0  second,  before  that 
all  four  satellites  are 
packed  together  in  a 
single  launch  vehicle. 
The  orbital  period  of 
the  final  orbits  is 


Fig.4.6  Distance  Between  Mother  and  First  Daughter  Satellites  10981.31  seconds,  so  the 


position  of  10981  s  is 

identical  with  the  point  of  0  second,  but  the  satellites  are  already  separated  500  km  apart. 
The  history  of  the  distance  between  the  mother  satellite  and  the  first  daughter  satellite  with 
respect  to  time  is  shown  in  Fig.4.6,  and  the  total  Av  required  for  the  four  satellites  is  shown 
in  Table  4.2. 
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4.4  STATIONKEEPING  MAINTENANCE 

The  analysis  to  this  point  is  based  on  the  Keplerian  orbit  theory12,  i.e.  there  is  no 
perturbation.  The  orbital  elements  are  very  sensitive  to  the  initial  condition.  Suppose  at  the 
time  t=l 098 1.31  sec.  the  mother  satellite  orbit’s  semi-major-axis  is  along  the  x  axis,  and  its 
semilatus  rectum  is  along  the  y  axis,  then  for  this  problem,  x  should  be  6978.137  km,  y=0, 
dx/dt=0,  dy/dt=8.77km/s;  for  the  first  daughter  satellite,  x=6962.313  km,  y=0.49977km, 
dx/dt=-0.5 1917km/s,  dy/dt=8.753km/s.  If,  somehow,  possibly  due  to  perturbations,  the  first 
daughter  satellite  has  a  small  initial  velocity  error  so  that  for  the  first  daughter  satellite 
dx/dt=-0.52842km/s,  dy/dt=8.787  km/s  then  the  distance  between  the  mother  satellite  and 
the  first  daughter  satellite  vs.  time  is  shown  in  Fig.4.7.  It  is  obvious  that  the  distance  can  not 
be  maintained  to  be  around  500  km.,  and  the  distance  is  diverging.  Even  though  the  first  and 

second  maneuvers  are 
performed  so  perfectly 
that  positions  and 
velocities  after  the 
second  maneuver  are 
exactly  as  required, 
with  perturbations 
(mainly  J2)  the  distance 
500  ±12  km.  can  not 
maintained  forever  (as 
indicated  in  Fig.4.4). 


Fig.4.7.  The  Distance  (km)  between  Adjacent  Satellites  vs 
Time  (sec.)  with  Small  Initial  Velocity  Error. 
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Fig.4.8  is  the  30  days 
simulation  result  of  the 
distance  between 

adjacent  satellites 

using  the  BG14.  From 
Fig.4.8,  it  is  seen  that 
with  perturbations  a 
secular  drift  appears  in 
the  transient  response; 
by  the  thirtieth  day  the 
distance  oscillates  around 

405±30  km.  Therefore,  some  kind  of  control  must  be  applied  to  the  formation  flight  system. 
There  are  many  papers  dealing  with  formation  control51,61,37,2;  most  of  them  are  approaches 
feeding  back  position  and  velocity  vector  errors.  In  this  paper  a  feedback  law  in  terms  of 
osculating  orbital  elements  is  studied.  By  comparing  the  situations  depicted  in  Fig.4.4  and 
Fig.4.8,  we  can  develop  a  control  strategy  so  that  the  control  effect  in  the  presence  of 
perturbations  (including  J2)  is  spent  principally  to  remove  the  secular  drift  of  the  separation 
distance  (Fig.4.8),  and  minimally  to  remove  the  amplitude  of  the  small  oscillations  (Fig.4.4). 
This  can  be  accomplished  with  the  feedback  of  the  variations  in  the  osculating  orbital 
elements  of  the  daughter  satellite  from  some  nominal  values.  With  this  philosophy,  control 
energy  can  be  conserved  as  compared  with  the  other  traditional  feedback  strategies,  e.g. 
feedback  of  position  and  velocity  error  components. 

Gauss’s  variational  equations  of  motion  provide  a  convenient  set  of  equations  relating 


Time(days) 

Fig.  4.8.  The  Distance  between  Adjacent  Satellites  with  Perturbations 
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the  effect  of  a  control  acceleration  vector  u  to  the  osculating  orbital  element  time 
derivatives52 


a  =  (2a2/h)  (  e  sin f  uf  +  (p/r)  uQ  )  (5a) 

e  =  [psin/  ur  +((p  +  r)cos/+re)  u0]//»  (56) 

i  =  [(r  cos0)  /  h ]  uh  (5  c) 

Q  =  [(r  sin0)/(6  sinl)]  uh  (5  d) 

co= — [-pcosfu  +(p+r)sinfuJ- - u.  (5e) 

he  6  sinl 

M=n  +[r\l(he)][(pcosf-2re)u  r~(p  +r)sin/M0]  (5 f) 


where  a  is  the  semi-major  axis;  e  is  the  eccentricity;  i  is  the  inclination;  Q  is  the  longitude 
of  the  ascending  node;  co  is  the  argument  of  the  perigee,  M  is  the  mean  anomaly.  We  define 
x=(a  e  t  Q  to  M)’  as  the  state  variable  vector  and  u=(ur  %  uj’  as  the  control  acceleration 
vector,  written  in  the  Local- Vertical-Local-Horizontal  frame,  Uj  points  radially  away  from 
the  Earth,  uh  is  aligned  with  the  orbit  angular  momentum  vector,  u9  is  orthogonal  to  both  u 
r  and  uh .  f  is  the  true  anomaly,  r  is  the  scalar  orbit  radius,  p  is  the  semilatus  rectum(see 
Eq.(4)),  0=G)+f,  h,  r|  and  the  mean  angular  velocity  n  are 

h  -\jp  y  T|  =\J\~e2  n  = ^  y/a 3 


Incorporating  the  J2  influence,  Eq.(5)  can  be  written  as  52 

x  =B(x)  u  +D(x)  (6) 


where 


Z)(x)=[0,  0,  0,  y,(— )2ncosl,  — /,(— )2n(5cos2i-l),  «+—/,(— )2qn(3cos2i-l)]r  (7) 

2  p  4  p  4  p 
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and 


B(x) 


(2a  2e  sin  f)/h 
(p  sin  f)/h 
0 
0 

-( p  cos f)/(he) 


(2  a2p)/(hr) 

[(p  +r)cos f~re\/h 
0 
0 

[( p+r)sinf]/(he ) 


(rcos0)//» 

(r  sin  0)/(/j  sin  i) 

-(r  sin  0  cos  i  )/(h  sin  i) 


T|(p  cosf~2re)/(he)  _  ["!](/?  +r)sinf]/(he) 


(8) 


Substituting  the  orbital  parameters  used  in  this  paper  into  Eq.(7) 

D(x)=[0  0  0  -5.218x10'®  -1.982x  10  7  n  +  1.919xl0'7]r  (9) 


The  elements  of  D  are  very  small,  and  if  the  formation  flying  system  is  under  control  the 
orbital  parameters  should  not  drift  far  away  from  the  nominal  values;  thus  D(x)  should  be 
very  close  to  Eq.(9).  Therefore,  D(x)  is  treated  as  a  minor  disturbance  instead  of  part  of  the 
plant  matrix. 

From  Eq.(8),  it  is  seen  that  the  system  is  nonlinear  and  time  variant,  so  a  control  law 
based  on  a  Lyapunov  function  is  applied.  If  the  osculating  orbital  elements  of  the  mother 
satellite  are  x„  the  required  osculating  orbital  elements  of  the  first  daughter  satellite  are  x2, 
then 

Ax=x2-X,  i.e.  x^jq+Ax  (10) 


Assuming  that  the  actual  osculating  orbital  elements  of  the  first  daughter  satellite  are  x2d, 
then 

6x=xw-x2  i.e.  x2d=x2+bx  (11) 
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We  define  a  Lyapunov  function  as 


V=-(a+be'a')  ■  bxT  bx 
2 


where  a> 0,  6>0,  a>0, 


V>0 


(12) 


then 


v  = 


bcte  a‘bx  Tbx+(a+be  a‘)bx  Tbx 

2 


(13) 


where 


bx=x2d-x2=x2d~xl-Ax  see  £^s.(ll)  and  (10) 

bx=x2d~xl  note  that  Ax  does  not  vary  with  time. 

substituting  Eqs.(6),  (7)  and  (8)  into  the  above  equation,  and  noticing  that  there  is  no  control 
for  the  mother  satellite,  we  get 


bx=B(x)u  +  [D(x2d)-D(x2)]  (14) 

Since  in  Eq.(7),  J2  and  Re  for  mother  and  first  daughter  satellites  are  the  same,  p,  n,  i,  and  r| 
are  almost  the  same,  and  D(x)  itself  is  very  small,  therefore,  D(x2d)-D(x1)  in  Eq.(14)  can  be 
ignored.  If  we  select 

«  =  -PCB  tB)'1B  Tbx  (15) 

where  P  is  a  scalar  value  used  to  adjust  the  feedback  gain,  and  substituting  Eq.(15)  into 
Eq.(14),  we  get 

bx  =  -pB(BTB)-lBTbx  (16) 
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After  substituting  Eq.(16)  into  Eq.(13)  there  results. 


V  =  -^bae~a,SxTSx-fi(a+be~a‘)5xTB(BTBylBT6x 


(17) 


<&=B(BTB)‘'BT  is  symmetric  and  semi-positive  definite,  and  Rank(T>)  <  Rank(B),  here  B 
is  6x3, the  eigenvalues  of  0  are  A,>0,  A2>0,  ^>0,  0,  0,  0.  There  must  be  an  orthogonal 
matrix  so  that 


<J>=rAr  T 


A,  o  oooo 


o  X2  o  o  o  o 


where  A=Tt<&T= 


0  A3  0 


0  0  0  0 
0  0  0  0 
0  0  0  0 


0  0 

0  0 
0  0 
0  0 


(18) 


Eq.(17)  can  be  written  as 


V=  -  6xT  0.5bae~a‘  T  I  Tr  6x  -  bxT  fi(a+bea')TAT  T5x 
=  -5 xT  [T  0.5  bae'a,I  Tt  +  T  P(a+6e  a‘)AT  t]5x 
=  -  6 xtT  [Q.5bae~a,I  +  $(a+bea,)A]T  T6x 
=  -  bxT  T  S  Tt  6x 


where 


S=[0.56ae'“7  +  P(a+6e“')A] 


(20) 


i.e. 


\0.5bae-at+$(a+be'at)L 


O.Sbae  "af+P(a+6e  _<W)A2 


0.5bae-at+f>(a+be'at)X.  0 


0  0.56ae  * 


Since  a>0,  b>0,  a>0,  (3>0,  and  A.,>0,  \2>0,  A3>0,  therefore  S  is  positive  definite,  and  so  is 

TETt.  If  8x  *0,  then  5xT  TETt  8x  >0, 


Mother 


Daughter 


p  (Bt  By'B1 


V=  -  6xT  T  a  Tt  dx  <  0 


That  means  the  control  law  described  by 
Eq.(15)  can  make  the  formation  flight 
system  asymptotically  stable.  It  is  obvious 
that  the  above  analysis  is  also  suitable  for 


the  control  of  the  distance  between  the 


Fig.4.9  Control  System  Diagram 


daughter  satellites,  e.g.  for  the  control  of 
the  second  daughter  satellite  to  maintain  the 


distance  between  it  and  the  first  daughter  satellite,  just  replace  x„  the  osculating  orbital 
elements  of  the  mother  satellite,  by  x  2,  the  requiredosculating  orbital  elements  of  the  first 
daughter  satellite;  replace  x2by  x  3  the  required  osculating  orbital  elements  of  the  second 
daughter  satellite,  and  replace  x2d  by  x  3d  the  actual  osculating  orbital  elements  of  the 
second  daughter  satellite,  in  Eqs.(10)  and  (1 1).  The  diagram  of  the  control  system  is  shown 


in  Fig.4.9. 
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Table  4.3.  Initial  Conditions  for  Mother  and  First  Daughter  Satellites 


Initial 

Condition 

Mother 

Satellite 

Daughter  Sat.l's 
Required  Elements 

Daughter  Sat.l’s 
Actual  Elements 

5x 

a  km 

•  10678.137 

10678.137 

10677.137 

-i 

e 

0.34650239 

0.34650239 

0.346534844 

0.000032453 

i  deg. 

83 

83 

82.9985 

-0.0015 

Q  deg. 

10 

10 

10 

0 

to  deg. 

10 

11.37 

11.37 

0 

M  deg. 

0 

0+39.5><n 

39.5xn 

0 

x  104 


Fig.4. 10  Transient  Responses  of  the  Semi-major 
Axis  and  the  Eccentricity 


4.5  NUMERICAL  SIMULATIONS 

With  various  perturbations, 
including  atmosphere  drag,  solar 
pressure,  the  Earth’s  magnetic  field, 
perturbations  from  the  Moon  and  other 
planets,  and  Jn’s  effect,  where  n=2  is 
the  most  important  contribution,  the 
osculating  orbital  elements  do  not 
remain  constant.  The  values  of  a,  the 
angle  between  the  semi-major  axes  of 
the  adjacent  satellites,  and  the  Pp,  the 
phase  distance  (see  Review  section)  at 
perigee  should  be  modified  when  a  or 
e  is  changed.  The  initial  conditions  for 


the  mother  and  first  daughter  satellites 
are  listed  in  Table  4.3.  When  P=2xl0'4  (see  Eq.(15)),  the  transient  responses  of  the  semi- 
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major  axis  and  the  eccentricity  for  ten  orbits  are  shown  in  Fig.4.10;  the  difference  of  the 
osculating  orbital  elements  between  the  mother  and  the  first  daughter  satellites  as  well  as  the 
control  efforts  of  the  first  daughter  satellite  are  shown  in  Fig.4.1 1.  From  these  figures  it  is 
seen  that  the  osculating  orbital  elements  converge  smoothly;  the  maximum  control  efforts 
(see  Fig.4.1 1,  g~I)  are  less  than  10'4  m/s2,  e.g.  if  the  mass  of  the  satellite  is  100  kg.,  the 
maximum  control  is  less  than  0.01  newton  (it  is  quite  small).  For  a  long  lifetime  satellite 
formation  flight,  this  should  be  reasonable.  The  main  purpose  of  this  study  is  to  maintain  the 
distance  between  adjacent  satellites;  the  distance  response  with  initial  error  and  various 


Fig.4.1 1 .  Difference  of  Orbital  Elements  Between  Adjacent  Satellites,  and  Control  Efforts 
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perturbations  for  30  orbits  is  shown  in 
Fig.4.12.  By  comparing  the  situations 
depicted  in  Fig.4.4  and  Fig.4.12,  we  can 
see  that  in  Fig.4.12,  the  initial  error  of  the 
distance  of  the  daughter  satellite  is 
gradually  diminishing,  after  several  orbits, 
the  distance  is  oscillating  around  500±12 


Fig.4. 1 2  Distance  Response  for  30  Orbits  km.,  just  like  the  situation  in  Fig.4.4,  which 
is  ideal  stationkeeping,  without  initial  error  or  perturbation.  From  Fig.4.12  it  is  also  seen 


that  although  the  control  effort  is 
small,  it  can  perfectly  compensate  the 
perturbations  (mainly  J2  effect).  The 
distribution  of  the  adjacent  satellites 
at  various  times  for  the  eleventh  orbit 
is  shown  in  Fig.4. 13. 

4.6  CONCLUSIONS 


x  107  m 


x  107m- 

Fig.4. 13  Distribution  of  Adjacent  Sat. 


A  strategy  for  maintaining  separation  distance  between  satellites  in  a  coplanar 
elliptical  orbiting  constellation  is  developed.  This  chapter  studies  the  implementation  of  this 
strategy.  This  strategy  would  be  implemented  by  two  maneuvers  that  would  cause  a  small 
angular  shift  in  the  directions  of  the  semi-major  axes  with  respect  to  the  semi-major  axis 
direction  of  a  “mother”  or  reference  satellite  which  is  also  included  within  the  constellation, 
and  a  certain  phase  difference  between  adjacent  satellites.  With  this  approach  for  Keplerian 
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type  orbits  and  nominal  alignment  of  the  satellites  along  the  orbit  track,  the  separation 
distance  between  adjacent  satellites  remains  within  a  few  percent  of  the  nominal  separation 
distance.  The  maneuvers  are  performed  at  the  perigee  point  in  the  orbit;  for  a  representative 
strawman  configuration  force-impulse  requirements  are  within  the  limits  for  proposed  pulse 
plasma  thrusters1.  It  is  clear  that  the  over-all  propulsion  requirements  would  be  far  less  than 
for  the  original  situation  in  a  Keplerian  orbit  without  the  semi-major  axis  angular  shift 
maneuver.  The  initial  separation  strategy  given  in  this  paper  results  in  the  minimum 
maneuver  energy,  since  the  work  done  by  the  thrusters  is  all  used  to  augment  the  orbital 
energy. 

For  stationkeeping,  a  control  law  based  on  the  feedback  of  the  tracking  errors  in  the 
osculating  orbital  elements  instead  of  feedback  of  the  traditional  Cartesian  position  and 
velocity  errors  is  presented.  One  of  the  benefits  of  this  feedback  law  is  the  removal  of  only 
the  secular  drift  caused  by  the  perturbations.  Therefore,  this  kind  of  control  can  save  energy; 
another  benefit  is  that  the  orbital  elements  which  do  not  have  tracking  errors  are  kept 
relatively  close  to  the  desired  values  during  the  maneuver. 

Since  the  model  is  nonlinear  and  time-variant,  the  control  law  is  based  on  a  Lyapunov 
function.  A  unique  Lyapunov  function  is  given  in  this  paper,  based  on  which  a  control  law 
is  established,  and  the  asymptotic  stability  of  the  formation  flight  system  is  strictly  proven. 
The  control  law  is  quite  simple,  so  it  is  easy  to  implement  in  engineering  practice. 

The  preliminary  analysis  reported  here  could  be  extended  to  more  complex 
geometrical  configurations  such  as  a  triangular  shaped  constellation  confined  to  the  orbit 
plane,  a  double  pyramid  configuration  with  three  spacecraft  in  the  common  bases  of  the 
pyramids  nominally  in  the  orbit  plane  and  two  spacecraft  at  the  vertex  points  of  the  pyramids 
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out  of  the  orbit  plane  of  the  base-plane  spacecraft.  Another  configuration  currently  proposed 
would  involve  four  in-plane  spacecraft  each  at  a  vertex  point  of  a  trapezoid.  For  the  more 
complex  in-plane  formations  both  along  track  and  radial  force  (thrusters)  would  be  required; 
for  the  out-of-plane  double  pyramid,  normal  forces  would  have  to  be  applied  to  the 
spacecraft  at  the  tips  of  the  pyramids. 
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V.  GENERAL  CONCLUSIONS  AND  RECOMMENDATIONS 


A  strategy  for  maintaining  separation  distance  between  satellites  in  a  coplanar 
elliptical  orbiting  constellation  has  been  developed.  This  strategy  can  be  implemented  by 
force-impulse  maneuvers  that  would  cause  a  small  angular  shift  in  the  direction  of  the  semi¬ 
major  axes  with  respect  to  the  semi-major  axis  direction  of  a  “mother”  or  reference  satellite 
which  is  also  included  within  the  constellation.  With  this  approach  for  Keplerian  type  orbits 
and  nominal  alignment  of  the  satellites  along  the  orbit  track,  the  separation  distance  between 
adjacent  satellites  remains  within  a  few  percent  of  the  nominal  separation  distance.  If  the 
maneuvers  are  performed  at  the  perigee  and/or  apogee  points  in  the  orbit  it  is  seen  for  a 
representative  strawman  configuration  that  force-impulse  requirements  are  within  the  limits 
for  proposed  pulse  plasma  thrusters. 

In  the  presence  of  perturbations  mainly  attributed  to  the  first  order  effects  of  the 
Earth’s  oblateness  for  the  highly  elliptical  strawman  configuration  orbit  the  results  are 
critically  dependent  on  the  amplitude  and  the  numerical  accuracy  in  calculating  the  semi¬ 
major  axis  shift  angle  for  a  nominal  along  track  separation  distance.  Without  subsequent 
corrections  a  secular  drift  is  observed  in  the  time  history  of  the  separation  distance  and  as  the 
time  increases  collision  or  near  collision  situations  can  exist. 

Additional  feedback  type  of  correctional  control  is  recommended  to  prevent  secular 
drifts  above  a  certain  level.  Two  types  of  stationkeeping  feedback  control  techniques  are 
considered  here:  (1)  an  application  of  linear  quadratic  regulator  (LQR)  theory  based  on  errors 
in  position  and  a  piecewise  adaptive  application  of  the  Tschauner-Hempel  equations  of 
motion  as  developed  for  elliptical  orbits,  and  (2)  also  based  on  a  Lyapunov  function  using 
osculating  orbital  elements.  For  the  Lyapunov  approach  the  asymptotic  stability  of  the 
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closed-loop  system  model  is  strictly  proven.  In  order  to  preserve  stationkeeping  fuel,  the 
feedback  type  of  control  should  be  initiated  only  when  the  errors  reach  a  certain  threshold. 

A  preliminary  deployment  strategy  is  introduced  based  on  near  Hohmann-type  of 
transfer  orbits.  Deployment  of  the  strawman  constellation  configuration  can  be  achieved  in 
one  orbit.  This  technique  was  selected  because  all  the  work  done  by  the  thrusters  is  used  to 
augment  the  energy  which  is  required  from  the  transient  orbit  to  the  final  orbit,  resulting  in 
near  minimum  maneuver  energy. 

The  preliminary  analysis  reported  here  could  be  extended  to  more  complex 
geometrical  configurations  such  as  a  triangular  or  trapezoidal  shaped  constellation  nominally 
located  within  the  orbital  plane,  or  a  double  pyramid  configuration  with  the  three  spacecraft 
in  the  common  bases  of  the  pyramids  nominally  in  the  orbital  plane  and  the  two  spacecraft 
at  the  vertex  points  of  the  pyramids  out  of  the  orbit  plane  of  the  base-plane  spacecraft. 
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VI.  IMPLICATIONS  FOR  FURTHER  RESEARCH 


1 .  Further  Parametric  Studies 

The  results  of  the  first  year  will  be  reviewed  and  extensions  implemented,  in  particular  the 
completion  of  deployment  maneuvers  based  on  the  solution  to  the  nonlinear  two  point 
boundary  value  problem  (TPBVP)  following  Pontryagin’s  principle.  This  work  is  in  progress 
at  the  end  of  the  first  year,  but  final  numerical  results  have  not  been  reported  in  this  volume. 
The  parametric  trade-off  studies  will  be  extended  to  study  the  effectiveness  of  both 
aerodynamic  drag  forces  and  solar  radiation  pressure  forces  for  formation  control  and  also, 
possibly,  for  attitude  control.  A  key  parameter  in  this  study  will  be  the  ballistic  coefficient 
of  each  satellite  in  the  constellation,  which  could  be  deliberately  altered  by  the  deployment 
of  extendible  vanes  from  some  or  all  of  the  satellites.  With  such  systems  the  use  of  hybrid 
control  strategies,  including  combinations  of  active  (thrusting)  control  as  well  as  the  semi¬ 
passive  control  afforded  by  the  orientation  of  the  vanes,  can  be  simulated  and  evaluated. 

2.  Implementation  of  Potential  Methods  to  Maintain  the  Formation. 

A  study  can  then  be  initiated  to  compare  full  active  propulsion  techniques  with  active 
propulsion  techniques  assisted  by  aerodynamic  /  solar  radiation  pressure  induced  forces. 
Among  full  active  propulsion  candidates  axe:  Pulse  plasma  thrusters  (PPT’s),  vaporizing 
ammonia,  resistojets  -  based  on  ammonia  or  hydrazine,  and  micro  -  propulsion  systems. 
Aerodynamically  /  solar  pressure  assisted  systems  could  consist  of  rotating  (movable)  drag 
panels  for  cross  orbit  track  and  along  the  track  or  a  combination  of  PPT’s  for  along  the  track 
and  aerodynamic  /  solar  pressure  panels  for  cross  orbit  track. 
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V 


3.  Use  of  Aerodynamics  /  Solar  Pressure  with  a  Deployable  Tether 

The  concept  of  using  rotating  aerodynamic  drag  /  solar  pressure  panels  pulling  against  the 
tether  to  maintain  cross  track  positioning  could  be  initiated.  HU  group  has  many  years 
experience  with  the  analysis  of  orbiting  tethered  systems  (Phase  A  and  B  of  Tethered  Shuttle 
Sub-satellite  System  for  NASA  and  Ball  Research,  followed  by  several  years  of  studying  the 
proposed  orbiting  tethered  antenna  /  reflector  system  for  AFOSR).  A  number  of  conference 
and  journal  publications  have  been  prepared  by  HU.  The  introduction  of  the  dynamics 
associated  with  a  deploying  or  retrieving  tether  into  the  constellation  system  dynamics  would 
involve  additional  modeling  and  simulation  requirements  which  would  need  to  be  added  to 
the  existing  software.  Initial  tether  modeling  could  be  based  on  a  relatively  short  tether 
assumed  to  be  massless  with  a  concentrated  end  mass  representing  the  panel.  Later  on 
modifications  could  account  for  the  effects  of  distributed  mass  along  the  tether,  a  much  more 
complex  modeling  problem. 

4.  Continuation  of  the  Use  of  Aerodynamic  Control  with  a  Deployable  Tether. 

Task  3  could  be  continued  to  include  more  complex  tether  modeling  as  required.  Different 
parametric  studies  would  be  conducted. 

5.  Analysis  of  Fixed  Formation  Flying  Based  on  A  Solid  Boom. 

The  solid  boom  could  represent  an  ultra  violet  hardened  tether,  a  miniature  Fairchild  type 
boom,  an  ultra  thin  “scissors  type”  boom,  or  a  thin  “tape  measure”  type  boom.  The  use  of 
such  a  boom  would  also  be  accompanied  by  solar  radiation  induced  forces  and  torques  in  the 
presence  of  associated  thermal  deflections,  which  may  be  important  depending  on  the  type 
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of  boom,  boom  thickness,  and  material  and  cross  sectional  characteristics. 

6.  Analysis  of  Spinning  Formation  Flying  Configurations. 

The  configuration  could  involve  subsatellites  as  end  bodies  connected  by  a  flexible  boom  or 
tether.  In  the  1970's  rotating  cable  connected  two-end  mass  configurations  were  studied  with 
the  aim  of  providing  a  certain  level  of  gravity  for  long  duration  manned  spaceflight 
applications.  This  concept  can  be  revisited  as  a  means  of  stabilizing  a  spinning  cluster 
satellite  configuration.  In  addition  to  the  study  of  the  complex  system  dynamics  (especially 
if  the  dynamics  of  the  end  bodies  are  considered),  a  study  of  the  effect  of  using  different 
materials  such  as  dental  floss  type  material,  Kevlar,  or  stainless  steel  should  also  be  included. 
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